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Abstract 

The  goal  of  this  research  is  to  build  up  a  logic  to  catch  and  track  incoming 
ASAT  weapons  by  using  space-based  onboard  optical  sensors.  The  satellite  orbit 
and  ASAT  trajectory  of  the  Chinese  test  were  generated  to  relate  the  research  to 
the  real  world  application.  These  position  and  velocity  values  are  used  to  generate 
simulated  observation  data  for  a  hypothetical  sensor  on  the  targeted  satellite.  These 
observation  values  are  assumed  to  be  true,  and  some  representative  amounts  of  error 
was  added  to  these  data.  Only  two  body  dynamics  are  considered;  drag  effect  and 
other  perturbations  are  neglected.  The  modern  orbit  determination  process,  least 
squares  method,  and  Monte  Carlo  techniques  are  used  to  calculate  the  estimated  orbit 
of  the  ASAT.  Standard  deviations  of  the  relative  position  of  the  ASAT  with  respect 
to  the  targeted  satellite  at  the  time  of  impact  are  calculated  for  different  sensors  with 
different  accuracy  and  data  collection  intervals. 
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Determining  the  capability  requirements 

FOR  A  SPACE-BASED  OPTICAL  SENSOR 
TO  DETERMINE  THE  TRAJECTORY  OF  AN 
INCOMING  ANTISATELLITE  WEAPON 

I.  Introduction 

1.1  Problem  Statement 

On  January  11,  2007,  the  People’s  Republic  of  China  (PRC)  conducted  its  first 
successful  antisatellite  (ASAT)  weapons  test  and  destroyed  its  own  Fengyun-lC,  a 
defunct  weather  satellite,  in  space.  FY-1C  was  launched  in  1999  and  was  orbiting 
the  Earth  in  a  polar,  low  Earth  orbit  (LEO)  at  an  altitude  of  about  537  miles  (865km), 
with  a  mass  of  about  750  kilograms.  China  used  a  modified  DON  FENG-21  road- 
mobile  Intermediate  Range  Ballistic  Missile  (IRBM)  with  a  600  kg  payload  as  an 
ASAT,  pictured  in  Figure  l.l-(b).  The  ASAT  missile  was  launched  from  China’s 
Xichang  Space  Center,  in  Sichuan  province,  shown  in  Figure  1.1- (a).  The  Chinese 
government  publicly  confirmed  that  they  had  conducted  the  test  on  January  23,  2007. 
Figure  l.l-(c)  represents  the  debris  cloud  created  just  after  impact.  The  calculations 
and  simulation  of  the  test  were  done  by  Dr.  T.  S.  Kelso  and  published  on  his  official 
website.  [6] 

The  test  had  two  major  and  concerning  results,  affecting  all  countries.  The  first 
of  which  was  the  potentially  damaging  space  debris  as  a  result  of  the  collision.  After 
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(a)  Xichang  Space  Center. 


(c)  Demonstration  of  Chinese  ASAT  test  by  AGI. 


Figure  1.1:  Chinese  ASAT  test.  [6] 
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the  test,  the  U.S.  Space  Surveillance  Network  tracked  and  cataloged  1337  pieces  larger 
than  10  centimeters  and  an  estimated  minimum  of  35000  additional  smaller  pieces 
were  declared.  The  resulting  debris  spread  from  200  kilometers  to  4000  kilometers 
and  began  to  endanger  almost  all  LEO  satellites  and  also  the  International  Space 
Station  (ISS),  orbiting  at  an  altitude  of  approximately  400  km.  The  Air  Force  Space 
Command  declared  that  the  space  debris  increased  the  collision  risk  for  about  700 
spacecraft.  According  to  the  optimistic  calculations  of  Dr.  Kelso  from  CclesTrack, 
after  100  years,  only  15%  of  the  total  debris  will  be  expected  to  have  decayed.  [6] 

On  the  other  hand,  and  more  relevant  to  this  thesis  research,  the  test  raised 
international  concerns  about  China’s  capability  and  intention  to  attack  satellites. 
This  was  the  first  destruction  of  a  satellite  with  an  ASAT  after  a  long  break  since  the 
Cold  War.  During  the  Cold  War  both  the  United  States  and  the  Soviet  Union  had 
conducted  such  ASAT  tests.  But  since  the  1980’s  neither  of  these  countries,  nor  any 
other  country,  has  intentionally  destroyed  satellites  in  space. 

In  such  an  arena  where  satellites  can  be  threatened  by  kinetic  energy  weapons, 
precautionary  actions  become  necessary  again.  Since  it  is  an  expensive,  time  con¬ 
suming  and  time-sensitive  process  to  put  a  payload  into  space  and  then  to  sustain  it 
in  order  to  get  the  advantage  of  the  ultimate  high  ground,  possessors  of  space  systems 
should  also  take  actions  to  protect  the  payloads.  Finally,  the  problem  appears  to  be 
seeing  and  avoiding  a  kinetic  kill  vehicle  intended  to  damage  or  destroy  the  satellite. 
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1.2  Method  of  Solution 


Optical  sensors  are  commonly  used  on  satellites  especially  for  attitude  determi¬ 
nation  and  control  subsystems  of  the  payloads.  Can  we  find  a  way  or  build  up  a 
logic  to  catch  and  track  the  incoming  ASAT  weapons  by  using  these  already  onboard 
optical  sensors?  The  goal  of  this  research  is  to  find  the  answer  to  this  question.  The 
satellite  orbit  and  ASAT  trajectory  of  the  Chinese  test  were  generated  to  relate  the 
research  to  the  real  world  application.  These  position  and  velocity  values  are  used  to 
generate  simulated  observation  data  for  a  hypothetical  sensor  on  the  targeted  satel¬ 
lite.  These  observation  values  are  assumed  to  be  true,  and  representative  amounts 
of  error  is  added  to  these  data.  A  modern  orbit  determination  process,  least  squares 
method,  and  Monte  Carlo  techniques  are  used  to  calculate  the  estimated  orbit  of  the 
ASAT.  Standard  deviations  of  the  relative  position  of  the  ASAT  with  respect  to 
the  targeted  satellite  at  the  time  of  impact  are  calculated  for  different  sensors  with 
different  accuracy  and  data  collection  intervals. 

1 . 3  Organization 

•  Chapter  1:  Introduces  the  problem  and  research  goals. 

•  Chapter  2:  Provides  a  comprehensive  problem  background  and  examines  what 
has  been  done  to  address  the  problem  area. 

•  Chapter  3:  Explicitly  details  the  problem  solving  approach.  Includes  an  exper¬ 
imental  method  description,  the  tools  and  techniques  developed,  and  approach 
verification  /  validation. 
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•  Chapter  4:  Discusses  and  provides  the  results.  Interprets  what  the  results  mean 
and  how  they  correlate  to  each  other. 

•  Chapter  5:  Covers  the  research  conclusions  and  major  result  trends.  Provides 
recommendations  for  future  research. 

•  Appendices:  MAT  LAB  code  supporting  the  simulations/experiments. 
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II.  Background  and  Literature  Review 


While  inspecting  the  background  of  the  research,  a  synopsis  of  the  previous  and 
current  implementations  on  space  weapons,  space  surveillance,  sensors  for  ballistic 
missile  defense  and  previous  research  will  be  beneficial. 

2.1  A  SAT  History 

An  antisatellite  weapon  can  be  defined  as  any  weapon  system,  whether  land- 
based,  sea-based,  airborne,  or  space-based,  which  is  specifically  designed  and  intended 
to  destroy,  damage  or  interfere  with  the  normal  functioning  of  space  objects.  Beside 
the  psychological  effects  on  the  enemy,  space’s  initial  military  use  was  reconnaissance. 
During  the  oversensitive  years  of  cold-war,  Sputnik  I  represented  the  idea  of  Russian 
superiority  in  space  technology  when  it  was  first  launched  in  1957.  The  launch  of 
Sputnik  also  triggered  the  desire  for  the  development  of  an  ASAT  weapon.  The 
first  official  project  of  the  Advanced  Research  Planning  Agency  (ARPA)  was  named 
Project  Defender,  covering  defense  from  both  satellites  and  ballistic  missiles.  The 
US  Air  Force’s  Air  Research  and  Development  Command  (ARDC)  signed  a  contract 
with  ARPA  for  the  “study  of  weapon  systems  to  combat  hostile  satellites”  in  1958. 
Later  NASA  was  involved  in  the  act  of  researching  ASAT  and  ballistic  missile  defense 
(BMD)  technology.  In  the  1960’s,  anti-satellite  capabilities  were  developed  as  part 
of  the  Soviet  space  defense  program.  They  began  development  of  a  limited  missile 
defense  of  Moscow,  which  employs  nuclear-tipped  interceptors.  [8] 


6 


The  very  first  ASATs  were  ballistic  missile  launched  weapons  with  either  non¬ 
nuclear  or  nuclear  warheads.  These  were  either  direct  hit-to-kill  devices,  or  satellites 
with  proximity  fuses  which  would  explode  using  debris  particles  to  destroy  the  target 
satellite.  In  the  case  of  using  an  ASAT  with  a  nuclear  warhead,  the  thermal  blast,  x- 
rays,  other  radiation  effects,  or  electro-magnetic  effects  would  be  the  kill  mechanism. 
An  early  Russian  ASAT  was  a  multi-staged  rocket  with  a  small  ground-controlled 
satellite  with  direct  hit-to-kill  capability.  A  self-guided  homing  vehicle  was  tried  using 
infrared  homing,  but  the  system  failed  several  times  in  testing,  was  not  successful  in 
the  60’s,  and  was  dropped.  [8] 

The  US  High  Altitude  Nuclear  Test  Program  was  a  study  to  determine  the 
effects  of  a  nuclear  warhead  explosion  in  space.  The  purpose  of  Project  Argus  was 
to  study  the  behavior  of  free  electrons  in  the  earth’s  magnetic  field.  The  US  military 
was  also  exploring  the  effect  of  nuclear  explosions  on  the  Explorer  IV  satellite,  which 
would  be  used  to  monitor  the  tests.  Three  nuclear  weapons  carried  by  rockets  were 
detonated  in  1958.  Project  Argus  showed  that  a  nuclear  explosion  in  space  generates 
high  energy  radiation  including  particles  from  the  explosion,  the  high  energy  electrons 
generate  radio  noise,  and  radio  transmissions  are  affected.  Also,  it  was  observed 
that  the  electrons  striking  the  metal  surfaces  of  satellites  can  damage  electronics. 
During  a  following  Fishbowl  project,  a  1.4  megaton  warhead  exploded  at  an  altitude 
of  248  miles  and  caused  considerable  interruption  with  Pacific  communications,  and 
destruction  in  power  systems  in  Hawaii,  and  damaged  three  satellites  in  orbit.  The 
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idea  of  using  nuclear  warheads  in  space  started  to  die  down  because  of  these  adverse 
effects  on  friendly  hardware.  [8],  [12] 

The  main  objective  of  Project  Bold  Orion  (rockets  launched  from  a  B-47  bomber, 
in  1959)  was  researching  the  feasibility  of  using  an  air-launched  ballistic  missile,  but 
the  project  was  extended  to  be  a  possible  ASAT  system.  Its  final  test  version  proved 
the  concept  and  approached  within  four  miles  of  the  Explorer  VI  satellite.  The  US 
Navy’s  Early  Spring  Project  was  proposed  to  mount  a  modified  Sparrow  air  to  air 
missile  on  a  sub-launched  ballistic  missile  (SLBM)  which  would  climb  to  the  target’s 
orbital  altitude,  wait  until  the  target  entered  the  engagement  zone,  and  then  destroy 
it  by  means  of  a  proximity  fuse  and  a  conventional  warhead.  Another  Navy  project 
Skipper  intended  to  launch  a  modified  Scout  rocket  from  a  ship  or  submarine  as  a 
kinetic  kill  ASAT  weapon.  The  Skipper  was  different  from  most  of  the  other  projects 
because  it  would  not  use  a  nuclear  warhead.  Neither  the  Early  Spring  nor  the  Skipper 
project  was  able  to  come  into  development.  [16] 

The  US  Army,  conducted  tests  using  the  Nike-Zeus  missile  system  which  was 
originally  developed  as  part  of  an  Anti-Ballistic  Missile  (ABM)  system.  The  first 
successful  US  anti-satellite  intercept  took  place  on  May  23,  1963  from  Kwajalein 
Island  in  the  Pacific  Ocean.  Throughout  the  duration  of  Project  Mudflap,  at  least 
eight  of  the  Nike  Zeus  ground-launched  missiles  were  fired  until  1966.  The  US  Air 
Force  deployed  Thor  rockets  which  were  modified  for  the  anti-satellite  mission  through 
Operation  Dominic  and  which  had  a  capability  of  carrying  a  1.5  megaton  yield  nuclear 
warhead  to  a  target  up  to  200  nautical  miles  high.  The  Dominic  project  conducted 


(a)  (b) 


Figure  2.1:  Figure  a,  b,  shows  Thor  and  Nike-Zeus  missiles 
respectively  .  [3],  [19] 

a  series  of  high  altitude  nuclear  tests  in  1962  and  later.  The  Thor  system  was  tested 
at  least  16  times  from  1964  to  1970,  prior  to  its  retirement  in  1976.  Figure  2.1  shows 
the  Nike-Zeus  and  Thor  missiles  used  as  ASAT  weapons.  [12],  [19] 

The  concept  of  using  Light  Amplification  by  Stimulated  Emission  of  Radiation 
(LASER)  and  Microwave  Amplification  by  Stimulated  Emission  of  Radiation  tech¬ 
nology  (MASER)  in  an  ASAT  role  of  attacking  the  satellite’s  sensors  and  electronics 
goes  back  to  the  early  1960’s.  The  United  States  was  aware  of  the  Soviet’s  improv¬ 
ing  development  of  particle  beam  ASAT  capability,  but  emphasis  remained  on  rocket 
powered  interceptors.  [12] 

Russia’s  main  ASAT  system  was  the  Co-Orbital  ASAT  system.  The  operation 
was  based  on  a  missile  armed  with  conventional  explosives.  The  1400  kg  ASAT 
interceptor  was  planned  to  be  placed  into  a  orbit  close  to  the  target  satellite’s  orbit 
and  then  it  maneuvered  to  get  close  to  the  target  within  one  or  two  orbits.  When 
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close  enough,  the  ASAT,  guided  by  its  onboard  radar,  was  directed  toward  the  target, 
exploding  and  destroying  the  target  using  fragmentation  effects.  The  project’s  initial 
testing  phase  was  between  1963  and  1972.  During  this  period  approximately  20 
launches,  seven  interceptions,  and  five  detonations  were  executed.  The  initial  tests 
confirmed  that  the  system  was  operational  from  altitudes  of  230  to  1,000  kilometers, 
and  the  system  was  declared  operational.  After  signing  the  Anti-Ballistic  Missile 
Treaty  in  1972,  the  Soviets  temporarily  ceased  tests,  but  in  1976  testing  of  the  Co- 
Orbital  system  resumed.  They  have  extended  the  operational  altitude  range  between 
160  km  and  1600  km,  minimized  the  attack  time  to  a  single  orbit  and  started  testing 
optical  and  infrared  sensors  in  addition  to  onboard  radar.  Testing  of  the  Soviet  Co- 
Orbital  ASAT  weapon  continued  from  1978  to  1982,  with  approximately  one  intercept 
per  year.  Although  it  has  not  been  tested  for  many  years,  the  system  is  thought  to 
be  currently  operational.  [8] 

The  Air  Force  ASAT  program,  Miniature  Homing  Vehicle  (MHV),  was  first 
mentioned  in  the  magazine  “Aviation  Week  &  Space  Technology”  in  March  1975. 
The  MHV  was  a  kinetic  energy  weapon  launched  from  an  F-15  and  guided  to  its 
target  by  an  infrared  sensor  mounted  in  its  nose.  This  weapon  was  composed  of  a 
small  two  stage  rocket  and  a  Miniature  Homing  Vehicle  (MHV)  which  would  destroy 
its  target  by  direct  impact  at  high  speed.  Launching  from  the  F-15  had  the  advantage 
of  being  able  to  bring  MHV  under  the  ground  track  of  its  target,  as  opposed  to  a 
ground-based  system,  which  must  wait  for  a  target  to  overfly  its  launch  site.  After 
a  long  period  of  development,  the  MHV  was  finally  tested  against  a  defunct  Army 
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Figure  2.2:  Miniature  Homing  Vehicle  launched  from  F-15.  [12] 

communications  satellite  in  1985,  successfully  destroying  the  target.  The  picture  at 
Figure  2.2  was  taken  during  one  of  the  MHV  tests.  [12] 

After  high  definition  imagery  from  satellites  in  low  Earth  became  very  accessible, 
the  US  reinvigorated  its  effort  to  find  a  way  to  neutralize  hostile  satellites.  The  US 
Army’s  Kinetic  Energy  Anti-Satellite  (KE  AS  AT)  program  started  in  1989  and  was  to 
provide  the  United  States  with  the  capability  to  destroy  hostile  satellites.  The  Kinetic 
Kill  Vehicle  (KKV)  would  be  launched  by  rocket  booster  to  strike  and  destroy  hostile 
satellites.  An  AS  AT  site  would  be  located  in  the  western  United  States  or  Pacific 
Ocean  area.  The  KE  ASAT  would  be  launched  when  the  target  approached  the 
firing  zone.  Then  the  KKV  would  separate  from  the  rocket  booster  and  make  course 
corrections  enabling  it  to  strike  the  satellite  and  disable  it  with  its  unique  debris 
mitigation  device.  After  this  the  KKV  would  re-enter  the  atmosphere  and  burn 
up.  The  KKV  already  had  been  designed,  developed,  integrated,  and  ground-tested 
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(a)  (b) 

Figure  2.3:  Kinetic  Energy  ASAT.  [12] 


successfully  and  was  supposed  to  have  become  mission  capabile  by  2000.  Figure  2.3 
depicts  the  KE  ASAT.  [12] 

Both  kinetic-kill  and  laser  ASATs  have  relative  advantages  and  disadvantages. 
Kinetic-kill  vehicle  systems  provide  easily  verifiable  destruction  of  a  satellite  and  are 
independent  of  weather.  Ground-based  lasers,  while  susceptible  to  adverse  weather, 
generate  less  space  debris  and  also  allow  for  a  covert  satellite  strike.  The  US  directed- 
energy  ASAT  system  centered  on  the  Mid-Infrared  Advanced  Chemical  Laser  (MIR- 
ACL)  which  was  developed  largely  in  1989  and  1990.  The  MIRACL  was  the  Erst 
megawatt-class,  continuous  wave,  chemical  laser  built  in  the  world.  The  US  Air 
Force  conducted  a  test  of  MIRACL  in  1997.  The  laser  was  directed  toward  a  satel¬ 
lite  orbiting  420  km  above  the  Earth.  The  MIRACL  laser  apparently  had  technical 
difficulties,  but  the  results  of  the  test  were  amazing.  A  lower-power  (30-watt)  laser, 
was  used  for  alignment  of  the  system  and  tracking  of  the  satellite,  but  it  was  observed 
that  this  lower-power  laser  was  powerful  enough  to  blind  the  satellite  temporarily, 
although  it  could  not  destroy  the  sensor.  In  the  1990’s  the  Soviets  also  developed  an 
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anti-satellite  laser  system,  and  their  system  was  considered  as  a  threat  to  satellites 
and  ballistic  missiles.  [12] 

Directed  Energy  Weapons  include  laser,  high-power  radio  frequency,  and  particle 
beam  technologies.  A  particle  beam  weapon  proposes  to  accelerate  charged  particles 
to  very  high  velocities.  High-energy  particle  beams  will  produce  high  surface  temper¬ 
atures  that  can  burn  out  the  satellite  electronics,  produce  high  surface  currents  that 
will  disrupt  sensitive  electronics,  or  produce  ions,  electrically  charged  particles  that 
will  disrupt  satellite  electronics  via  radiation  effects.  On-the-ot her- hand,  it  would 
be  difficult  and  expensive  to  place  particle  beam  weapons  in  orbit.  Many  tons  of 
material  must  be  lifted  and  a  complex  device  must  be  constructed  under  free-fall 
conditions.  Electronic  signal  manipulation  is  another  major  class  of  ASAT  weapons 
effort.  The  signal  to  the  satellite  can  be  disrupted  with  a  very  high,  electronic- 
competing  signal.  Traditional  satellite  components  are  also  becoming  smaller  and 
lighter.  This  may  eventually  permit  the  launch  of  parasitic  microsatellites  which  can 
maneuver  close  enough  to  the  target  satellite  to  disrupt  or  destroy  it.  In  recent  days, 
the  US  has  begun  working  on  several  systems  including  the  Experimental  Spacecraft 
System  (XSS-11),  the  Near-Field  Infrared  Experiment  (NFIRE),  and  the  Space-based 
Interceptor  (SBI)  programs.  [12] 

China’s  ASAT  test  of  11  January  is  the  sort  of  capability  available  to  any  country 
which  has  IRBMs  or  satellite  launch  vehicles  and  a  long-range  radar  system,  such 
as  Japan,  India,  Pakistan,  Iran,  and  even  North  Korea.  Many  countries  now  use 
space  systems  for  military  and  intelligence  purposes.  In  addition  to  the  US  and 
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Russia,  for  example,  several  European  countries,  Israel,  India  and  Japan  also  maintain 
reconnaissance  satellites  in  LEO,  and  these  are  all  vulnerable  to  ASAT  missiles.  [12], 

[U], 

2.2  Space  Surveillance 

2.2.1  Overview.  For  the  most  part,  the  initial  intent  of  space  surveillance 
sensors  was  to  provide  warning  of  a  strategic  missile  attack.  However,  the  rising 
number  of  satellites  created  a  requirement  to  watch  all  satellites  during  their  lifecy¬ 
cle,  including  launch  and  decay,  in  order  not  to  confuse  them  with  hostile  missiles. 
Eventually,  the  space  surveillance  operations  started  to  separate  from  missile  defense 
operations  with  the  increase  of  the  military  importance  of  space. 

Both  optical  and  radar  systems  are  used  as  satellite  tracking  systems  and  they 
mostly  use  the  latest  and  most  expensive  sensor  technologies.  Most  optical  sensors 
are  dependent  on  reflected  sunlight  or  emitted  infrared  energy  to  track  a  satellite. 
On-the-other-hand,  active  optical  sensors,  illuminating  a  target  with  coherent  laser 
radiation  are  being  used  in  some  recent  applications.  By  illuminating  a  target  with 
laser  radiation,  these  systems  can  image  satellites  that  are  not  reflecting  sunlight.  Ac¬ 
tive  illumination  also  provides  measurement  of  the  range  to  the  target.  Ground-based 
radar  systems  have  been  used  since  the  late  1950s  for  space  surveillance  applications. 
Radars  have  the  advantage  of  being  able  to  track  the  target  any  time  and  uninter¬ 
ruptedly  during  cloudy  conditions  as  compared  to  optical  sensors.  Today,  using  the 
advanced  technology  of  large  phased  array  radars  (LPAR),  a  great  variety  of  opera- 
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tions  including  satellite  tracking,  missile  early  warning,  and  guiding  interceptors  as 


apart  of  ABM  systems  can  be  accomplished. 


2.2.2  United  States  Space  Surveillance  Network.  One  part  of  United  States 
Space  Command  (USSPACECOM)’s  mission  is  to  detect,  track,  catalog  and  identify 
man-made  objects  orbiting  Earth.  These  include  active  and  defunct  satellites,  rocket 
bodies,  and  debris.  As  described  by  the  USSPACECOM  itself,  space  surveillance 
accomplishes  the  following: 


•  “Predict  when  and  where  a  decaying  space  object  will  re-enter  the 
Earth’s  atmosphere; 

•  Prevent  a  returning  space  object,  which  to  radar  looks  like  a  missile, 
from  triggering  a  false  alarm  in  missile-attack  warning  sensors  of  the 
US  and  other  countries; 

•  Chart  the  present  position  of  space  objects  and  plot  their  anticipated 
orbital  paths; 

•  Detect  new  man-made  objects  in  space; 

•  Produce  a  running  catalog  of  man-made  space  objects; 

•  Determine  which  country  owns  a  re-entering  space  object; 

•  Inform  NASA  whether  or  not  objects  may  interfere  with  the  space 
shuttle  and  Russian  Mir  space  station  orbits.  [1]” 

The  Space  Surveillance  Network  (SSN)  accomplishes  these  efforts  via  US  Army, 
Navy  and  Air  Force  operated  ground-based  radars  and  optical  sensors  at  25  sites 
worldwide.  The  SSN  started  tracking  Sputnik  I  and  is  still  observing  and  tracking 
space  objects. 

The  SSN  has  tracked  more  than  24,500  space  objects  orbiting  Earth  since  the 
launch  of  Sputnik  I  in  1957.  Most  of  these  objects  have  decayed  entering  the  Earth’s 
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atmosphere,  and  the  SSN  tracks  more  than  8,000  orbiting  objects  now.  Beside  satel¬ 
lites,  rocket  body  pieces  of  10  pounds  and  space  objects  as  small  as  10  centimeters 
in  diameter  can  be  tracked  by  SSN.  The  number  of  operational  satellites  is  approxi¬ 
mately  seven  percent  of  the  total  number  of  8,000;  all  of  the  other  objects  are  debris. 
The  SSN  predicts  the  space  objects’  orbits  and  checks  the  objects  at  an  instant  rather 
then  continuously  following  them  because  of  the  limited  number  of  sensors  and  other 
limited  capabilities  of  the  network.  The  US  SSN  radar  sensors  and  their  held  of  view 
at  500  km  altitude  is  pictured  in  Figure  2.4. 


Figure  2.4:  Space  Surveillance  Network  Radar  Sensors  and 

their  held  of  view  at  500  km  Altitude.  [11] 


The  SSN  consists  of  the  following: 

Phased-array  radars  have  no  moving  mechanical  parts,  and  the  radar  energy  is 
directed  electronically.  Using  this  advantage,  these  radars  can  track  multiple 
satellites  simultaneously  and  scan  large  areas  of  space  in  a  very  short  time. 
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For  instance,  the  AN/FPS-85  phased-array  radar  at  Eglin  AFB  in  Florida  is 
composed  of  almost  6,000  transmitter  antennas  and  20,000  receiver  antennas 
and  can  track  objects  from  just  above  the  horizon  to  very  close  to  the  zenith  over 
an  azimuth  of  120  degrees.  It  can  track  space  objects  up  to  40,000  kilometers 
in  range. 

Conventional  radars  composed  of  tracking  and  immobile  antennas  which  operate 
in  bistatic  mode.  Bistatic  mode  means  one  antenna  transmits  a  pulse  and 
another  receives  the  reflected  signals.  The  Naval  Space  Surveillance  System 
(NAVSPASUR),  operating  with  conventional  radars,  is  a  network  of  three  trans¬ 
mitting  and  six  receiving  radars  providing  continuous  observation  and  detection 
of  space  objects  crossing  the  continental  United  States. 

Ground-Based  Electro-Optical  Deep  Space  Surveillance  System  (GEODSS) 
consists  of  three  operational  sites  at  Socorro,  New  Mexico;  Maui,  Hawaii;  and 
Diego  Garcia,  British  Indian  Ocean  Territories.  GEODSS  combines  the  func¬ 
tions  of  telescope,  low-light-level  television  and  computers  to  perform  the  surveil¬ 
lance  missions.  Each  site  has  three  telescopes.  An  example  complex  is  seen 
in  Figure  2.6.  The  system  only  operates  at  night  but  has  the  capability  to 
detect  objects  10,000  times  dimmer  than  the  human  eye  can  detect.  Since  it 
is  an  optical  system,  it  is  adversely  effected  by  cloud  cover  and  local  weather 
conditions.  GEODSS  can  track  objects  as  small  as  a  basketball  at  a  distance 
of  more  than  20,000  miles  in  space.  It  has  a  critical  role  of  tracking  deep  space 
objects  including  geostationary  satellites.  The  location  of  the  GEODSS  sites 
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with  optical  track  capability  and  their  coverage  at  500  kilometers  can  be  seen 


at  Figure  2.5. 


Figure  2.5:  Locations  of  GEODSS  Sensors  and  their  field  of 
view  at  500  km  Altitude.  [11] 


Figure  2.6:  GEODSS  Site.  (Reference:  space.kursknet.ru) 


All  these  different  types  of  sensors,  located  at  different  SSN  sites  such  as  Maui, 
Eglin,  Thule,  and  Diego  Garcia  collect  up  to  80,000  satellite  observations  each  day. 
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Each  station  transmits  its  data  directly  to  USSPACECOM’s  Space  Control  Center 
(SCC)  by  means  of  satellite,  ground  wire,  microwave,  and  phone.  The  SCC  in 
Cheyenne  Mountain  Air  Station  is  the  nucleus  of  the  SSN  where  the  enormous  amount 
of  data  flows.  The  Cheyenne  Mountain  command  center  is  seen  in  Figure  2.7.  The 
SCC  uses  computer  aided  systems  to  process  SSN  information  and  accomplish  the 
space  surveillance  and  space  control  missions.  The  NAVSPACECOM  is  the  site  for 
the  Alternate  SCC  (ASCC).  [11],  [17]. 


Figure  2.7:  Cheyenne  Mountain  Command  Center.  [13] 


2.2.3  Russian  and  European  Space  Surveillance  Network.  The  Russian  space 
surveillance  system  uses  an  the  early-warning  radar  network  and  is  operated  by  the 
space-surveillance  division  of  the  3rd  Army.  The  network  also  includes  the  Krona 
system  at  Zelenchukskaya  in  the  North  Caucasus  and  Nakhodka  on  the  Far  East. 
The  main  optical  observation  station  that  monitors  objects  on  high-altitude  orbits  is 
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called  Okno  and  is  located  in  Nurek,  Tajikistan.  Okno  can  detect  objects  at  altitudes 


of  up  to  40,000  km. 


Figure  2.8:  The  FGAN  Tracking  and  Imaging  Radar  (TIRA) 
at  Wachtberg,  Germany.  [14] 

Although  ESA  and  European  countries  with  space  monitoring  capabilities  are 
strongly  dependent  on  initial  object  and  orbit  information  provided  by  USSPACE- 
COM,  European  countries  have  the  capability  to  track  the  Earth  orbital  environment 
up  to  and  beyond  GEO  altitudes.  Main  coordination  of  the  systems  is  done  by  the 
European  Space  Agency.  The  ESA  Space  Debris  Telescope,  the  French  ROSACE/- 
TAROT  system,  and  the  UK  PIMS  sensors  can  detect  GEO  objects  well  below  the 
stated  USSPACECOM  size  threshold  of  1  m  in  diameter.  The  GRAVES  receiver  at 
Apt,  France;  the  FGAN  Tracking  &  Imaging  Radar  (TIRA)  at  Wachtberg,  Germany 
(shown  at  Figure  2.8);  Phased-array  surveillance  radar  and  tracking  radars  at  Fyling- 
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dales,  UK;  Norwegian  Globus  II  radar;  British  Radar  at  Fylingdales,  UK;  the  French 
GRAVES  system;  the  Chilbolton  radar  located  in  Winchester, UK  and  the  European 
Incoherent  Scatter  Radar  are  some  of  the  powerful  radars  that  are  used  for  space 
surveillance  and  early  warning  operations.  [14] 

2.3  Ballistic  Missile  Defense  Platforms  and  Sensors 

For  an  effective  ballistic  missile  defence,  a  group  of  sensors  which  are  a  combi¬ 
nation  of  land,  sea  or  space  based  sensors  should  be  incorporated  to  detect  and  track 
a  threat  missile  through  its  trajectory.  Only  a  worldwide  sensor  coverage  network 
can  track  a  missile  during  all  boost,  midcourse  and  terminal  phases. 

Defense  Support  Program  (DSP)  Satellites  are  orbiting  the  earth  approximately 
35,780  kilometers  over  the  equator  in  a  geosynchronous  orbit.  The  system 
provides  global  coverage  for  early  warning,  tracking  and  identification  using  in¬ 
frared  sensors  to  detect  heat  from  missile  and  booster  plumes  against  the  Earth’s 
background.  Recently,  in  addition  to  their  primary  mission  of  ballistic  missile 
defense,  DSP  satellites’  infrared  sensors  started  to  be  used  in  an  early  warning 
system  for  natural  disasters  like  volcanic  eruptions  and  forest  fires. 

Space  Based  Infrared  System  (SBIRS)  program  is  the  follow-on  capability  to 
the  Defense  Support  Program  (DSP).  DSP  has  has  been  used  for  more  than 
30  years  as  an  early  warning  systems  for  ballistic  missile  launches.  Now  the 
goal  of  the  US  is  to  provide  transition  from  DSP  to  SBIRS  without  any  gap  in 
the  ABM  defense  system.  The  SBIRS  program  currently  includes  two  Geosyn- 
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chronous  Earth  Orbit  (GEO)  satellites,  two  Highly  Elliptical  Orbit  (HEO)  pay- 
loads  installed  on  GEO  satellites,  and  ground  systems  deployed  around  the 
world.  SBIRS  provides  early  warning  of  ballistic  missile  attacks  and  high  pre¬ 
cision  information  for  ABM  systems  to  intercept  and  destroy  threat  missiles. 
The  system  is  currently  under  development  by  the  US  Air  Force  and  a  proposed 
future  constellation  is  pictured  in  Figure  2.9. 


Early  Warning  Radars  (EWR)  are  aiming  to  determine  the  final  destination  of 
threat  missiles  more  precisely.  Ground-based  radars  located  in  California, 
Alaska  and  overseas  are  being  upgraded  by  the  US  in  order  to  be  used  more 
effectively  with  the  developing  ballistic  missile  defense  system. 

X-Band  Radars  are  capable  of  searching,  detecting  and  tracking  missiles.  They 
can  also  distinguish  between  warheads  and  decoys.  One  of  the  latest  version  of 
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X-Band  radars  is  currently  being  constructed  by  the  US  Missile  Defense  Agency. 
It  is  an  X-band  radar  mounted  on  a  moveable  semi-submersible  oil  drilling 
platform.  Sea  based  X-band  radar,  shown  in  Figure  2.10,  will  be  home-ported 
in  Adak,  Alaska  and  will  increase  the  overall  success  of  ballistic  missile  defense 
systems  with  its  capability  to  move  throughout  ocean  areas  for  operations  and 
testing. 


Figure  2.10:  Sea  based  X-band  radar.  [1] 


THAAD  Radars  are  the  main  sensors  of  the  Terminal  High  Altitude  Area  Defense 
(THAAD)  system.  This  valuable  subsystem  of  the  BMD  system  is  rapidly 
transportable  and  forward-deployable.  The  upgraded  version  will  have  the 
capability  to  intercept  and  destroy  ballistic  missiles  inside  or  outside  the  atmo¬ 
sphere  while  they  are  in  their  final,  or  terminal  phase  of  flight. 

SPY-1  Radar  is  mounted  on  Aegis  cruisers  and  destroyers  and  a  part  of  the  BMD 
Agency’s  sea  based  ABM  system.  Aegis  Destroyers  with  S-Band  phased  array 
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SPY-1  radars  detect  and  track  Intercontinental  Ballistic  Missiles  and  report 
track  data  to  other  systems  and  sites  of  BMD  sysytem. 

Forward  Deployable  Radars  (FDR)  can  even  be  air  transportable  and  provide 
additional  sensor  capability  for  tracking  hostile  missiles.  Forward  Deployable 
Radars  are  intended  to  be  placed  at  sites  close  to  the  launch  area  of  ballistic  mis¬ 
siles  where  they  can  obtain  more  accurate  tracking  data  very  early  and  transfer 
this  data  to  friendly  interceptors  providing  additional  time  for  more  successful 
defensive  intercepts.  [1] 

2.4  Previous  Research 

The  Space-Based  Visible  Program  was  aiming  to  accomplish  the  technology 
demonstration  of  space  based  space  surveillance  with  the  Midcourse  Space  Experiment 
satellite.  The  satellite  was  launched  in  1996  into  a  near  sun-synchronous  LEO  orbit 
with  an  onboard  visible-band  electro-optical  camera  designed  at  Lincoln  Laboratory. 
In  1997,  technology  demonstration  was  successfully  achieved  by  gathering  optical 
information  on  various  space  objects.  Then,  the  space-based  visible  sensor  associated 
with  the  program  started  to  be  used  as  a  part  of  US  Space  Surveillance  Network  with 
a  proved  capability  at  least  as  accurate  as  GEODSS  sensors  of  the  network.  [7] 

In  a  paper,  C.  B.  Chang  examined  the  problem  of  ballistic  trajectory  estimation 
with  angle-only  measurements.  Earth  gravity  was  assumed  as  the  predominant  and 
only  force  on  the  long  range  ballistic  missile  trajectory.  An  iterative  least  squares 
algorithm  was  used  to  determine  the  state  of  the  trajectory  with  angle  only  measure- 
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merits.  Both  ground-based  and  space-based  sensors  were  considered.  The  estimation 
algorithm  described  in  this  paper  is  emphasized  as  applicable  to  state  estimation  prob¬ 
lems.  [5] 

In  the  paper  “Ballistic  Missile  Track  Initiation  From  Satellite  Observations”, 
written  by  Murali  Yeddanapudi,  Yaakov  Bar-Shalom,  Krishna  R.  Pattipati,  Somnath 
Deb,  an  algorithm  is  presented  to  track  a  ballistic  missile  in  the  initial  phase,  out 
of  atmosphere,  using  line-of-sight  measurements  from  one  or  more  moving  sensors, 
typically  mounted  on  satellites.  Results  of  the  estimation  problem  were  considered  as 
non-satisfactory  because  of  the  poor  target  motion  capability  when  using  the  Gauss- 
Newton  iterative  least  squares  minimization  algorithm  for  estimating  the  state  of 
a  nonlinear  deterministic  system  with  nonlinear  noisy  measurements.  The  major 
problem  with  this  approach  was  caused  by  strong  dependence  on  initial  guess.  A  more 
sophisticated  Levenberg-Marquardt  method  was  used  in  place  of  the  simpler  Gauss- 
Newton  algorithm  and  robust  new  methods  for  obtaining  the  initial  guess  in  both 
single  and  multiple  satellite  scenarios  were  developed.  The  sensor  was  considered  to 
be  both  passive  and  active.  Monte  Carlo  simulation  studies  on  some  typical  scenarios 
were  performed,  and  the  results  indicate  that  the  proposed  estimators  are  efficient. 
[15] 

Most  of  the  studies  about  ASAT  systems  and  anti-ASAT  systems  are  classified 
or  have  a  restricted  distribution.  As  an  example,  the  AFIT  thesis  titled  ‘Protection  of 
a  High  Valued  Space  Asset’  describes  a  current  gap  at  situational  space  awareness  and 
offers  an  onboard  sensor  system  as  a  preferred  solution.  Candidate  active  and  passive 
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sensors  are  revealed.  As  a  short  term  protection  solution  for  space  based  high-valued 
assets,  employing  optical  sensors  for  each  payload  is  recommended.  Restrictions  on 
this  research  prevent  its  exact  recommendations  from  being  published.  However,  just 
knowing  the  type  of  recommended  sensor  helped  to  focus  this  thesis.  [4] 

2.5  Chapter  Summary 

ASAT  systems  have  been  commonly  developed,  tested,  and  prepared  for  use 
since  1950’s.  The  technology  developed  rapidly  after  the  first  space  launches.  With 
the  recent  Chinese  test,  it  is  proven  one  more  time  that  KE  ASAT  weapons  can  eas¬ 
ily  be  used  to  threaten  LEO  satellites.  Improved  applications  of  rocket  technology, 
guidance  systems  and  more  effective  warheads  make  the  space-based  assets  more  vul¬ 
nerable.  On-the-other-hand,  space  surveillance  and  ballistic  missile  defense  efforts, 
briefly  described  through  the  chapter,  shows  that  precautions  against  ASAT  weapons 
should  be  composed  of  worldwide  network  of  systems,  including  highly  accurate  sen¬ 
sors.  Previous  researches  emphasize  that  even  the  US  has  some  gaps  in  defending 
their  space  assets  against  ASATs.  Also,  there  are  a  lot  of  countries  that  have  satellites 
but  don’t  have  space  surveillance  sensors.  It  is  obvious  that,  if  the  decision  is  based 
upon  the  information  coming  from  external  data  sources,  on  most  occasions  there  will 
not  be  enough  time  to  maneuver  the  satellite  defensively.  Finally,  the  idea  of  having 
a  sensor  onboard  the  space  asset  with  an  early  warning  capability  has  matured.  If 
the  satellite  can  watch  the  earth,  especially  the  suspicious  areas  and  detect  and  track 
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a  possible  incoming  ballistic-missile,  it  might  have  enough  time  to  defeat  the  missile 
and  survive. 
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III.  Methodology 


In  this  section,  the  methodology  of  estimating  the  ASAT  missile’s  trajectory 
based  on  a  space-based  optical  sensor’s  measurements  will  be  analyzed.  Initially,  the 
orbit  of  the  FY-1C  satellite  and  the  trajectory  of  the  ASAT  missile  were  created  using 
two  body  dynamics.  The  state  vectors  including  position  and  velocity  for  each  of 
them  were  created  throughout  the  flight  time  of  the  ASAT  missile  at  different  time 
steps.  These  data  assumed  to  be  true  and  true  observations  were  generated  using 
these  state  vectors.  Then,  some  noise  was  added  to  those  observations  and  simulated 
real  world  observations  were  computed.  The  least  squares  estimation  filter  which  was 
gathering  these  observation  data  as  input  and  generating  estimated  state  of  the  ASAT 
missile  at  epoch  time,  was  set  up.  Finally,  the  most  probable  converged  estimated 
state  at  epoch  time  was  propagated  in  time  and  the  position  of  the  ASAT  missile 
with  respect  to  the  target  at  the  time  of  impact  was  found.  The  flow  of  this  process 
is  explained  in  Figure  3.1. 


ADO  SOM:  NOISE 


Figure  3.1:  Analysis  Flow  Diagram. 


3.1  Assumptions 

•  In  the  Principia  Newton  formulated  the  law  of  gravity  beside  his  three  laws 


of  motion.  The  law  of  gravity  is  expressed  in  Equation  3.1,  where  G  is  the 


universal  gravitational  constant,  F  is  the  gravitational  force  on  mass  m  caused 
by  mass  M  and  r  is  the  vector  pointing  from  M  to  m.  Any  set  of  two  bodies 
attract  each  other  with  a  force  proportional  to  the  product  of  their  masses 
and  inversely  proportional  to  the  square  of  the  distance  between  them.  This 
principal  is  the  basic  law  for  the  motion  of  the  satellites  and  planets. 

Few, TV  =  -F^r-  (3.1) 

However,  when  the  relative  motion  of  two  bodies  in  real  world  is  considered,  it 
has  been  observed  and  proved  that  some  perturbations  are  involved  in  to  the 
dynamics.  A  perturbation  can  generally  be  described  as  a  deviation  from  the 
expected  behavior.  Some  forces  caused  by  atmospheric  drag  and  lift,  thrust, 
solar  radiation,  magnetic  and  relativistic  effects,  the  Earth’s  and  the  Moon’s 
non-perfect-spherical  shapes,  and  forces  due  to  additional  space  objects  all  cause 
deviations  during  the  relative  motion  of  the  two  bodies.  The  two-body  prob¬ 
lem  was  described  to  simplify  the  equation  of  motion  for  only  two  bodies  with 
assumptions  that  the  masses  of  the  bodies  can  be  concentrated  at  their  centers, 
and  there  is  not  any  force  other  than  gravitational  force  acting  on  the  system 
along  the  line  joining  the  centers  of  the  bodies.  To  accomplish  the  initial  study 
about  this  topic,  the  two-body  equation  of  motion  (EOM)  is  applied  to  the  mo¬ 
tion  of  the  satellite  and  the  ASAT  weapon  targeting  the  satellite.  The  other 
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perturbations  affecting  their  relative  motion  are  neglected  and  assumed  to  be 
easily  added  in  future  studies.  [18] 

•  The  simulated  Chinese  ASAT  weapon  is  assumed  to  be  a  modified  but  unguided 
ballistic  missile.  The  missile  does  not  make  any  correction  on  its  track  to  the 
targeted  satellite. 

•  The  optical  sensor  on  the  satellite  can  establish  line  of  sight  measurements  with 
an  accuracy  of  up  to  1/1000  Arcseconds. 

•  The  missile  can  be  tracked  throughout  its  entire  flight  path  including  launch 
and  impact.  The  atmospheric  affects  on  the  optical  observations  like  reflection, 
refraction  or  the  attenuation  of  the  light  are  neglected. 

•  Since  the  curvature  of  the  earth  is  not  considered  in  the  calculations,  it  is  as¬ 
sumed  that  the  line  of  sight  between  the  sensor  and  the  ASAT  missile  is  not 
obstructed  by  the  earth. 

•  As  described  in  the  two-body  problem,  neither  the  satellite’s  nor  the  ASAT 
missile’s  orbits  are  influenced  by  atmospheric  affects.  The  trajectory  of  the 
ASAT  missile  is  assumed  to  be  always  exo-atmospheric. 

•  It  is  supposed  that  the  satellite  is  able  to  calculate  its  position  and  velocity 
vectors  with  respect  to  the  geocentric  equatorial  coordinate  system  at  every 
thousandth  of  a  second  and  its  calculations  are  assumed  to  be  free  of  error. 
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3. 2  Scenario 


The  actual  Chinese  AS  AT  test  conducted  on  January  2007  is  represented  in  the 
scenario.  The  Chinese  defunct  weather  satellite  FENG  YUN-lC’s  orbital  parameters 
were  available  as  open  public  information.  Using  those  data  the  state  vector  of  the 
satellite  Xsat—  [x,y,z,  Vx,Vy,Vz}T ,  composed  of  its  position  and  the  velocity  vectors 
with  respect  to  the  inertial  frame,  were  created.  The  state  vector  of  the  satellite  at 
the  impact  time,  5:28  p.m.  EST  on  January  11,  2007,  was  also  determined.  Since  we 
have  the  information  about  the  location  of  the  launch  area,  which  was  the  Chinese 
Xichang  space  complex,  we  knew  the  start  and  the  end  points  of  the  ASAT  missile’s 
trajectory.  Applying  approximate  launch  parameters  for  an  IRBM  which  is  capable 
of  this  mission,  and  the  location  of  the  two  end-points,  the  trajectory  which  will  best 
fit  in  these  data  was  created.  In  general  there  would  be  two  possible  impact  points 
in  the  missile  trajectory  but  the  problem  was  set  up  such  that  the  impact  occurs 
at  the  missile’s  apogee.  In  these  calculations,  the  closest  interception  between  two 
orbital  tracks  is  selected  as  the  impact  point  and  the  state  vector  of  the  ASAT  missile 
Xref=  [x,  y,  z,  Vx,  Vy,  VZ]T  was  created.  The  state  vectors  for  both  the  satellite  and 
the  ASAT  missile  were  generated  for  different  time  intervals  throughout  the  entire 
missile  flight  time  of  approximately  8  minutes.  These  data  or  the  software  program 
used  to  create  these  data  are  not  demonstrated  in  this  thesis  because  they  are  straight¬ 
forward  and  simple  to  reproduce.  These  vectors  were  assumed  to  be  perfect  and  used 
to  calculate  the  simulated  observations  achieved  by  the  onboard  sensor. 
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The  scenario  is  such  that  we  assume  the  sensor  on  the  satellite  was  tracking  the 
earth  continuously  and  had  the  locations  of  the  possible  threat  launch  sites  stored 
in  its  database.  Whenever  it  caught  a  threat  it  would  begin  to  make  observations 
and  calculate  the  position  of  the  missile  at  each  observation  time  step.  After  having 
enough  data  it  would  calculate  the  estimated  track  of  the  ASAT  missile  and  compare 
it  with  its  own  orbit  in  order  to  estimate  a  possible  intersection.  The  number  of 
observations  required  to  make  an  estimate  depends  on  the  accuracy  of  the  sensor 
and  the  time  intervals  between  the  observations.  Once  the  sensor  obtained  enough 
observations  it  would  propagate  the  missile  trajectory  and  calculate  the  standard 
deviation  of  the  miss  distance  between  the  two  orbits  at  the  possible  intersection 
time.  Being  aware  of  this  possible  intersection  error  ellipsoid  would  provide  the 
satellite  the  chance  to  make  a  last  ditch  manoeuver  to  defeat  the  ASAT  missile. 

3.3  Dynamics 

As  mentioned  in  the  assumptions  section,  in  this  research  numerical  integra¬ 
tions  of  the  two-body  problem  were  used  to  calculate  the  state  vectors  of  the  targeted 
satellite  and  the  ASAT  missile.  Initially,  the  state  vectors,  composed  of  position  and 
velocity  vectors  with  respect  to  the  geocentric  equatorial  coordinate  system,  were  cre¬ 
ated  and  they  were  assumed  to  be  the  “ truth Then  simulated  observation  measure¬ 
ments  were  created  based  on  this  true  data  and  some  statistically  relevant  additional 
error.  A  brief  explanation  of  the  two-body  problem  is  considered  to  be  helpful. 
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The  two-body  problem  is  a  result  of  Newton’s  second  law,  and  two  assumptions 
that  (1)  the  bodies  can  be  represented  as  point  masses  and  (2)  there  is  no  other  force 
acting  on  them  other  than  the  gravitational  force.  To  be  able  to  describe  the  relative 
motion  of  the  two  bodies,  an  inertial  frame  should  be  described.  Newton  described 
this  inertial  frame  as  “in  its  own  nature,  without  relation  to  anything  external,  remains 
always  similar  and  immovable.  ”  There  is  an  ongoing  discussion  on  his  inertial  frame 
definition  but  it  can  be  assumed  as  almost  inertial  and  then  proceed  to  explain  the 
two-body  problem.  [18] 


z 


Figure  3.2:  Relative  Motion  of  Two  Bodies. 

In  the  book  named  Fundamentals  of  Astrodynamics  [18]  two  bodies  with  masses 
M  and  m  are  considered  and  pictured  as  in  Figure  3.2.  Their  position  is  defined  with 
respect  to  an  inertial  frame  {X' ,  Y' ,  Z')  and  a  secondary  non  rotating  reference  frame 
(. X ,  Y,  Z),  which  is  parallel  to  the  inertial  frame  and  originated  at  the  center  of  mass 
M.  The  vector  r  is  defined  as  r  =  rm  —  vm  ■ 
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Applying  Newton’s  second  law  in  the  inertial  frame  the  dynamics  equations 


become 


mr 


m 


MrM 


GMm  r 

/r>2  /T> 

GMm  r 


These  equations  are  simplified  and  written  as 

GM 

^  m 

Gm 

r  m  =  —  r 


(3.2) 

(3.3) 


(3.4) 

(3.5) 


Subtracting  Equation  (3.5)  from  Equation  (3.4)  the  vector  differential  equation  of  the 
relative  motion  for  the  two-body  problem  is  obtained  as 


G(M  +  m) 
r  = - — - r 


(3.6) 


When  the  situation  in  which  a  smaller  object  like  an  artificial  satellite  or  a  ballistic 
missile  orbiting  a  planet  is  considered,  the  mass  of  the  orbiting  m  will  be  very  small 
when  compared  to  the  central  body  M ,  so  it  can  be  assumed  that 


G(M  +  m )  «  GM 


(3.7) 


34 


Finally,  when  the  gravitational  parameter  /i  =  GM  is  substituted  in  Equation  (3.6), 
the  convenient  form  of  the  two-body  equation  of  motion  is  obtained  in  Equation  (3.8). 


l-i 

r  +  —  r  =  0 


(3.8) 


3.4  Filter 

3. 4-1  Least  Squares  Estimation.  In  the  two-body  problem,  six  classical  or¬ 
bital  elements  define  the  orbit.  There  should  be  at  least  six  observations  in  order  to 
calculate  these  six  orbital  elements.  Although  it  can  be  assumed  that  the  dynamics 
are  perfect,  the  observations  include  errors  and  in  most  applications  there  are  more 
observations  than  six.  A  German  mathematics  student  Carl  Freidrich  Gauss  dis¬ 
covered  the  Principal  of  Maximum  Likelihood  and  the  Least  Squares  method,  solved 
these  orbit  determination  problems  and  succeeded  to  determine  the  orbit  of  asteroid 
Ceres.  Ceres  was  first  observed  by  Piazzi  but  lost  after  a  few  observations.  Gauss 
calculated  the  orbit  that  passes  as  close  as  possible  to  all  observation  points  obtained 
by  Piazzi  and  the  lost  asteroid  Ceres  was  discovered  using  the  orbit  determined  with 
least  squares  estimation.  [20] 

3. 4-1.1  The  Principle  of  Maximum  Likelihood.  Initially,  it  can  be 
assumed  that  in  order  to  find  the  wanted  value  of  xo,  N  independent  measurements 
of  Xi  were  obtained  using  different  equipment,  and  each  data  had  standard  deviation 
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(Tj.  The  joint  probability  can  be  calculated  as  the  product  of  the  individual  Gaussian 


distributions. 
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(3.9) 


Since  the  true  value  Xo  in  Equation  (3.9)  is  always  unobtainable,  an  estimate  x  of  the 
wanted  true  xq  value  is  replaced  for  it. 
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(3.10) 


The  value  of  x  which  is  closest  to  the  truth  can  be  obtained  by  maximizing  the 
probability  in  Equation  (3.10)  with  respect  to  x.  And  the  maximum  of  that  equation 
can  be  obtained  when  the  argument  of  the  exponential  is  minimized. 


d  (xi  —  x)2 

tic  2a? 

2=1  L 


(3.11) 


This  step  of  the  principal  is  alternately  named  as  the  method  of  least  squares.  [20] 


3. f.  1.2  Linearization  of  Dynamics  .  The  equations  of  motion  used 
in  this  thesis  are  especially  based  on  the  two-body  dynamics.  State  vectors  a;  of  a 
satellite  and  an  ASAT  missile  were  created  and  estimated.  The  state  vectors  of  the 
dynamic  system  were  composed  of  position  and  velocity  vectors.  As  Wiesel  describes 
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in  his  text,  the  equations  of  motion  determine  the  evaluation  of  the  state  vector  with 


time  as  the  differential  equations 


di=eM  (3-12) 

or  the  closed-form  solution  can  be  written  as  relating  the  actual  state  to  the  initial 
state  and  time. 


x(t)  =  h(x(t0),t )  (3.13) 

Both  of  these  equations  can  be  used  to  calculate  the  changes  of  the  true  state 
xq,  the  estimated  state  x,  and  the  nearby  trajectories.  It  can  be  assumed  that  the 
estimate  of  the  true  trajectory  will  be  very  close  to  the  true  trajectory  and  express 
their  relation  as  x  =  x0  +  Sx.  Substituting  this  expression  into  the  equations  of 
motion  (3.12),  expanding  g  in  a  Taylor’s  series,  and  ignoring  second  and  higher  order 
terms  the  equation  becomes 

^  =  g(x0  +  Sx,t) 

~  g(x0,t) +  Vxg(x0,t)5x  +  0(2)  (3.14) 

Using  Equation  (3.14)  the  Equations  of  Variation  are  expressed  as 
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—  Sx  =  Ait )  Sx 
dt  v  ' 


(3.15) 


where 

Mt)  =  v,  g,o(t)  (3.16) 

Recalling  the  fact  that  equations  of  variations  are  linear  ordinary  differential  equa¬ 
tions,  Wiesel  derived  and  defined  the  State  Transition  Matrix  <f>  as 

8x(t)  =  Q(t,to)  8x(t o)  (3-17) 

where  <f>  satisfies  the  equation 

^  $(Mo)  =  A(t)  $(t,t0)  (3.18) 

with  initial  conditions 

$(t0,t0)  =  I  (3.19) 

In  Equation  (3.19)  the  term  /  represents  the  identity  matrix.  In  linear  systems 
<f>  is  used  to  propagate  the  state  in  time;  however,  in  nonlinear  systems  $  matrix 
propagates  the  small  variation  of  the  state  5x  in  time.  Whenever  the  equations  of 
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motion  are  available  in  closed  form,  the  state  transition  matrix  for  each  time  step  can 
be  obtained  by  numerical  integration.  [20] 

3-4- 1-3  Linear  Least  Squares.  Linear  least  squares  method  is  used 
to  estimate  the  state  of  a  linear  dynamical  system  at  an  epoch  time  t0  with  using 
the  observations  Zi(ti)  taken  at  time  t,t.  Each  observation  vector  Zi(ti )  is  assumed 
to  be  independent,  unbiased  and  has  an  associated  instrumental  covariance  Qi.  If 
observation  data  is  linearly  related  to  the  system  state  at  the  measurement  time, 
observation  relation  can  be  expressed  as 


Zi(ti)  =  Hi  x(ti)  +  et 


(3.20) 


in  which  e%  is  the  true  error  of  the  observation.  If  we  insert  the  system  dynamics  into 
this  relation  and  substitute  the  definition  T*  =  Hi  $(£j,t0)  we  get  a  more  compact 
observation  relation. 


Ziiti)  =  H.l  ${ti,t0)x(tQ)  +  e< 

Zi(tj)  =  Tix(t0)  +  e*  (3.21) 


Following  the  method  in  Wiescl’s  text,  the  total  data  vector  z,  the  observation 
matrix  T  and  the  instrumental  covariance  matrix  Q  can  be  created 
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Zi 

Z2 

Zn 


(3.22) 


T\ 

t0) 

t2 

= 

to) 

tn 

HNQ(tNi  to) 

Qi  0  •••  0 
0  Q2  ■■■  0 


0  0  •  ■  •  Qv 


(3.23) 


(3.24) 


for  N  number  of  observations.  And  after  some  derivations  Wiesel  defines  the  estimate 
of  the  state  vector  x(t0)  and  the  covariance  matrix  Px(t0)  at  an  epoch  time.  [20] 


x(to)  -  (TtQ-1T)-'TtQ-1z 

(3.25) 

=  (TrQ-lT)-' 

(3.26) 
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34.1.4  Nonlinear  Least  Squares.  Nonlinear  least  squares  method  is 
much  more  appropriate  when  solving  the  real  world  problems  compared  to  the  linear 
least  squares,  because  in  real  world  most  of  the  dynamics  and  the  observation  relations 
are  nonlinear.  Also,  the  solution  of  the  two-body  problem,  used  in  this  thesis,  is  highly 
nonlinear.  Linearization  of  these  nonlinear  system  dynamics  can  be  established  with 
some  assumptions.  Again,  following  Wiesel’s  method,  assuming  system  dynamics  are 
available,  they  can  be  expressed  in  two  forms  as 


dx 

dt 


(3.27) 


x{t)  =  h(x(t0),t) 


(3.28) 


where  Equation  (3.27)  represents  the  case  when  we  have  the  equations  of  motion  and 
Equation  (3.28)  is  the  explicit  solution  of  the  dynamics  as  a  function  of  time  and  the 
initial  conditions.  Assuming  that  the  dynamics  are  deterministic,  linearization  of  the 
dynamics  about  a  reference  trajectory  xref 


Sx(t )  =  <£>(£,  t0)  dx(t0)  (3.29) 

should  be  valid.  The  reference  trajectory  xref  is  expected  to  be  close  to  the  estimated 
trajectory  x  which  is  used  instead  of  the  unobtainable  true  trajectory  xq,  and  5x  is 
expected  to  be  small.  The  5x  amount  of  changes  will  correct  the  reference  trajectory 
xref  and  turn  it  into  the  closest  estimate  x  of  the  true  trajectory  x0. 
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The  observations  will  be  nonlinearly  related  to  the  current  state  at  that  time. 
This  can  be  expressed  as 

Zi(ti)  =  G(x(ti),ti)  (3.30) 

where  z(ti)  are  the  measurements  taken  during  the  observation  at  time  tt.  In  real 
world  applications  all  measurements  include  some  error.  If  we  were  able  to  find  a 
perfect  instrument,  it  would  observe  the  true  state  Xq  and  generate  the  true  observa¬ 
tion  data  zq.  The  vector  of  actual  measurements  with  true  error  is  defined  as  z  and 
the  imperfect  observed  state  is  called  x.  Assuming  that  true  error  in  the  data  goes 
to  zero  when  the  error  in  the  state  approaches  to  zero,  and  also  assuming  that  they 
are  both  small  enough,  the  true  error  in  the  actual  data  becomes 


Z  —  Zq  =  G(x,t )  —  G(xo,t) 

(3.31) 

G(xo  +  Sx,t)—  G(xo,t ) 

(3.32) 

8GC  ,  , 

(3.33) 

where  x  =  xq  +  5x,  e  ~  r  (in  which  r  represents  the  residuals )  and  the  matrix  in 
the  Equation  (3.33)  relates  the  errors  in  the  state  to  the  reference  trajectory.  The 
residuals  can  be  represented  as 

n  =  Zi-  G(xref(ti),ti)  (3.34) 
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Similar  to  the  linear  case,  H  matrix  is  for  the  nonlinear  case  is  defined  as 


SG 

Hi  -fix  ) 


(3.35) 


Using  the  fact  that  these  residuals  are  linearly  related  to  Sx,  which  changes  the  ref¬ 
erence  trajectory  to  into  an  estimated  trajectory,  and  Sx  propagates  as  shown  in 
Equation  (3.29),  the  residuals  equation  becomes 


HiSx(ti )  =  Hi^(ti,t0)Sx(t0) 

(3.36) 

TiSx(t0 ) 

(3.37) 

Finally,  the  solution  for  the  nonlinear  least  squares  estimation  can  be  written  as 


8x(t0)  =  (TTg-1r)-1rTg-1r 

(3.38) 

pSx  =  (TTg_1r)_1 

(3.39) 

where  Sx  is  the  small  variances  to  the  reference  trajectory  and  P$x  is  the  covariance 
matrix.  The  general  form  of  the  covariance  matrix  in  Equation  3.39  is  a  symmetric 
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matrix 


P 1 1  Pi-2  •  •  •  Pin 
Pl2  P22  '  '  '  P2N 

P  = 

Pin  P2N  ■  ■  ■  Pnn 

where  the  diagonal  components  Pa  are  the  variances  representing  the  a2  values  and 
the  off  diagonal  components  represent  the  covariances.  In  a  special  case  where  all 
the  measurements  are  statistically  independent  of  each  other  the  covariance  matrix 


(3.40) 


becomes  [20] 

a2  0  •  •  •  0 
0  a  1  ■  ■  ■  0 

P  =  (3.41) 

0  0  •  •  •  a% 

and  as  Wiesel  explains  [20],  the  estimated  trajectory  is  determined  as 


x(t0)  =  xref(to)  +  5x(t0)  (3.42) 

3-4-2  Observation  Geometry.  The  purpose  of  this  thesis  is  to  analyze 
whether  only  one  optical  sensor  designated  for  self  defense  can  estimate  a  threat’s 
trajectory  based  only  on  its  observed  measurements.  Although  we  don’t  have  any 
real  world  data  to  be  used  as  the  observations,  measurements  including  simulated 
errors  were  generated  using  the  assumed  true  trajectories.  As  mentioned  previously, 
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the  orbit  of  the  satellite  and  the  trajectory  of  the  ASAT  missile  were  created  via 
numerical  integration  of  the  equations  of  motion  throughout  the  flight  time  of  the 
missile.  The  state  vectors  composed  of  position  (x,  y,  z )  and  velocity  (14,  Vy,  Vz)  vec¬ 
tors  for  each  time  interval  4  were  generated.  The  initial  conditions  were  selected  to 
fit  in  the  scenario  and  achieve  the  impact  at  the  given  time  using  the  closest  inter¬ 
section  point.  In  order  to  do  the  integration  in  MAT  LAB,  the  equations  of  motion 
were  written  in  matrix  form. 


X  =  B(X)X 


(3.43) 
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(3.44) 


The  optical  sensor  on  the  satellite  is  assumed  to  be  able  to  achieve  a  line  of  sight 
measurement  to  an  incoming  ASAT  missile.  The  exact  position  of  the  satellite  rsat 
is  gathered  from  the  central  computer  of  the  payload  and  it  assumed  to  be  true  for 
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the  observation  time  fj.  Given  that  the  sensor  is  aware  of  the  possible  AS  AT  missile 
launch  locations  and  the  approximate  launch  parameters  of  a  possible  incoming  missile 
we  can  say  that  the  sensor  can  calculate  the  approximate  position  of  the  ASAT  missile 

T asat • 


FY-1C 


1 

Figure  3.3:  Observation  Geometry. 

The  observation  to  be  used  is  defined  as  the  cosine  of  the  angle  between  the 
position  vectors  of  the  satellite  and  the  ASAT  missile.  If  a  is  the  observation  angle, 
z  =  cos  a  is  the  observation  data  measured  by  the  sensor.  This  angle  is  defined  as 
a  theoretical  angle  to  be  used.  Some  other  angles  like  the  angle  between  the  lines 
connecting  the  center  of  the  earth,  the  satellite  and  the  ASAT  could  also  be  used. 
The  observation  is  not  perfect  and  the  observed  angle  a  is  related  to  the  true  angle 


46 


a  between  the  true  position  vectors  as  a  =  a  +  <70,  and  represented  in  Figure  3.3. 
Then,  the  observation  can  be  written  as 

z  =  cos  a  =  cos  (a  +  Sift)  (3.45) 

where  cos  a  =  Zo  is  the  true  observation  without  error  and  SiJ)  can  be  defined  as 
the  la  accuracy  of  the  sensor;  later  this  accuracy  will  be  multiplied  by  a  zero  mean, 
discrete-time  white  Gaussian  random  number  to  generate  error  based  on  the  fact  that 
real  instruments  typically  exhibit  such  statistics  [9].  The  error  in  the  measurement 
can  be  expressed  as 

\8z\  =  |  cos  a  —  cosd|  =  |  cos  (a  +  Sift)  —  cos  o:|  (3.46) 

and  expanding  the  cos  (d  +  Sip)  term  the  observation  error  becomes 

\8z\  =  |  cos  a  cos  Sip  —  sin  a  sin  8ip  ■—  cos  a \  (3.47) 

Since  the  angle  8ip  is  assumed  to  be  small  enough  and  with  small  angle  assumptions 

cos  S'lfi  ~  1 
sin  S'i/j  ~  8ip 


Equation  (3.47)  becomes 
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(3.48) 


\8z\  =  |  —  5'0sino;| 

Squaring  both  sides  and  substituting  the  zo  =  cos  a  we  get 

(h,)2  =  (<^)2  (1  -  4)  (3.49) 

Finally,  applying  estimation  operator  as  described  by  Wiesel  [20]  the  instrumental 
variance  for  each  observation  at  time  tt  can  be  obtained 

(a2)2  =  M2  (1  -  ^)  =  Q  (3.50) 

in  which  ay,  is  the  variance  of  the  sensor  as  given  by  the  manufacturer  or  determined 
by  experiment. 

In  this  analysis  initially,  zq  =  cos  a  is  calculated 


zo 


cos  a  = 


F sat  '  ^asat 


|  f sat 1 1  f  asat \ 

•E sat-E asat.  4“  V  satV  asat  4~  Zsa^^asaf 

'\/^'sat  4”  V Sat  4~  ^saty/^asat  4"~  Vasat  4"~  Zasaf 


(3.51) 


and  some  error  5z,  as  defined  in  the  Equation  (3.49),  is  added  to  each  true  observation 
in  order  to  get  the  simulated  real  world  measurements.  In  order  to  calculate  erroneous 
simulated  real  world  measurements  z%  at  each  observation  time  tu  the  S ij}  term  in 
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Figure  3.4:  Gaussian  Distribution  of  the  Error.  [2] 


Equation  (3.49)  is  multiplied  with  a  random  number,  generated  in  MAT  LAB  using 


Gaussian  distribution.  Observation  data  for  the  entire  flight  time  of  the  missile  is 


obtained  as 


Zi 

ZOi 

ci\j0-  -4) 

-2 

= 

-02 

+ 

c2y/(l-4) 

Zn 

Z0N 

cN^{l-zlN) 

(3.52) 


where  N  is  the  total  number  of  measurements  and  ci,  C2,  ...Cat  are  the  random  numbers 
created  based  on  Gaussian  distribution  with  a  mean  of  S'ljj,  which  is  shown  in  Figure 
3.4.  [2] 


3-4-3  Filter  Processing.  The  objective  of  this  research  is  to  determine  the 
relative  position  of  an  unguided  ballistic  missile,  modified  as  an  ASAT  weapon,  with 
respect  to  the  targeted  satellite  at  the  time  of  impact.  This  problem  can  be  solved 
with  a  combination  of  estimator,  dynamics  model  and  data  observations.  The  solution 
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of  this  problem  is  attempted  to  be  achieved  using  a  space-based  onboard  optical 
sensor,  two-body  dynamics,  and  least  squares  estimation.  The  nonlinear  dynamical 
model  to  estimate  the  state  of  ASAT  in  this  problem  can  be  shown  to  be 


X(ti)  =  h(X(t0),ti )  (3.53) 

where  X(t0)  is  the  state  of  the  system  at  the  epoch  time  t0  and  h  is  an  analytic 
solution  to  the  equations  of  motion.  The  observation  relation  wit  the  state  is  defined 
in 

z(ti)  =  G[X(ti)]  +  v(ti)  (3.54) 

where  z{ti)  is  the  observation  made  at  time  ti  and  v(ti)  is  an  independent  zero-mean 
discrete-time  white  Gaussian  noise  with  a  variance  of  Q.  Further,  in  our  case 


=  G(Xref)  =  COS  £  =  XsatXref  +  VsatVJf  +  Z™iZre^_  (355) 

\/ Xlat  +  Vsat  +  Zsat  J Xref  +  Vref  +  Zref 


and  as  defined  in  the  nonlinear  least  squares  section  Hi  is 


Hi  =  —  (Xref(ti),ti)  (3.56) 


If  we  define  rsati  =  y/xjat.  +  y'jati  +  z2sat.  and  rrefi  =  yjx2ref.  +  y2refi  +  z%f.,  the  H 
matrix  for  each  time  tt  becomes 
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(3.57) 
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(3.58) 


In  accordance  with  the  algorithm  for  nonlinear  least  squares,  as  mentioned  by 
Wiesel  [20],  the  state  vector  should  be  propagated  to  the  observation  time  £j.  In 
order  to  solve  the  state  transition  matrix  d>  simultaneously  with  the  state  X  their 
expression  in  the  two-body  equation  of  motion  must  be  rearranged  and  written  in 
matrix  form  to  be  numerically  integrated  in  MAT  LAB.  Recalling  the  two-body 
equation  of  motion 

r  +  ^  r  =  0  (3.59) 

and  the  composition  of  the  ASAT’s  state  vector  with  respect  to  the  geocentric  equa¬ 
torial  coordinate  frame 
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and  combining  them  in  matrix  form  we  get 


X  =  B(X)X 


(3.61) 
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(3.62) 


Additionally,  to  get  the  state  transition  matrix  for  the  observation  time  ti:  the 
equations  of  variation  can  be  numerically  integrated  as 


^  $(Mo)  =  A(t)$(t,t0) 


(3.63) 


where  A  —  V/ 


dh 

dh 

dfi 

dh 

dh 

dh 

dx 

dy 

dz 

dVx 

avy 

avz 

dh 

dh 

dh 

dh 

dh 

dh 

dx 

dy 

dz 

d\ 4 

dVy 

dVz 

dh 

dfs 

dh 

dh 

dh 

dh 

A(t)  = 

dx 

dy 

dz 

dVx 

9Vy 

avz 

dh 

9f± 

dh 

dfi 

dh 

dfi 

dx 

dy 

dz 

dVx 

dVy 

avz 

dh 

dh 

<9/5 

dh 

dh 

dh 

dx 

dy 

dz 

dVx 

dVy 

avz 

dfe 

dfe 

dfe 

dh 

dfe 

dh 

dx 

dy 

dz 

dVx 

dVy 

avz 

Substituting  the  the  equations  we  get 
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(3.65) 


and  when  we  get  the  partial  derivatives,  we  have 
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(3.66) 
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where  /  is  the  identity  matrix  and  <p  represents  the  null  matrix.  As  in  Wiesel’s 


algorithm  Arr  can  be  created  as 


A 


rr 


/i  _|_  3 fix2  3/ixy 


3  fixz 


3/ixy 

y-5 


I  3  fiy2 
y  3  ^*5 


3  nyz 


3/ixz 


3  yyz 


(3.67) 


Now,  if  we  write  the  equations  of  variation  in  matrix  form  we  get 


- 

- 

— 

“ 

$11 

$12  ‘ 

'  $16 

Al 

A12  ' 

7^16 

$11 

$12  ' 

'  $16 

$21 

$22  ' 

'  $26 

= 

^21 

A22  ' 

A26 

$21 

$22  ' 

'  $26 

$61 

$62  ‘ 

1 

CO 

.  co 

^61 

A§2  ' 

1 

CO 

CO 

$61 

$62  ' 

- 1 

CD 

CO 

(3.68) 


In  order  to  numerically  integrate  Equation  (3.68)  using  MAT  LAB  it  has  to  be  mod¬ 
ified.  Also,  the  state  vector  and  the  state  transition  matrix  should  be  integrated 
simultaneously.  Finally,  the  matrix  form  that  can  be  integrated  is  set  up  as 
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The  time  interval  of  the  numerical  integration  can  be  any  time  period  through¬ 
out  the  flight  time  of  the  ASAT  missile  with  selectable  time  steps.  The  initial 
conditions  should  be  IqX6  for  the  state  transition  matrix  <f>  and  any  selected  state 
on  the  reference  trajectory  for  the  initial  state  of  the  ASAT  missile.  Returning  to 
Wiesel’s  algorithm,  the  filter  calculates  residuals  as  rt  =  Zi  —  G(X).  Hi  and  the 
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observation  matrix  Tt  =  HiQ  are  calculated  for  each  time  step  of  the  observations. 
The  covariance  of  the  correction  P$x  and  the  correction  vector  5X(to)  at  epoch  time 
are  calculated  by  the  filter  using  running  sums  at  each  iteration. 


P&x 

5X(t0) 


PsxY,T?Q7lr- 


(3.72) 

(3.73) 


Given  in  Equation  3.72  and  Equation  3.73,  y^lTfQi  1Tij  should  be  invertible  to 
hncl  SX. Finally,  the  filter  corrects  the  reference  trajectory  state  at  epoch  time. 


-l 


Xref+l(^o)  Xref(to)  +  SX{to)  (3.74) 

After  Ending  the  corrected  state  at  epoch  time  the  filter  uses  this  new  state  as 
the  initial  condition  of  the  ASAT  trajectory.  The  filter  repeats  the  iteration  until 
it  reaches  the  convergence  criteria.  The  convergence  criteria  can  come  from  the 
covariance  matrix.  It  is  determined  that  if  the  correction  to  each  element  of  the  state 
vector  is  smaller  than  half  of  the  square  root  of  the  corresponding  diagonal  element  of 
the  covariance  matrix  there  is  no  need  to  continue  iterating.  To  explain  convergence 
criteria  symbolically,  for  the  components  of  the  last  computed  state  vector  if  all  of 
the  following  are  simultaneously  true 
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8x 

< 

0.5^/ '  P$x 

(U) 

5y 

< 

0.5^  Psx 

(2.2) 

Sz 

< 

0.5^  Psx 

(3,3) 

5VX 

< 

0.5a/  Psx 

(4,4) 

6Vy 

< 

0.5^  Psx 

(5,5) 

svz 

< 

0.5  y/P^c 

(6,6) 

then  the  filter  stops  iterating,  takes  the  last  computed  state  at  epoch  time  to  as  initial 
conditions,  propagates  the  state  in  time  and  finds  the  position  and  velocity  of  the 
AS  AT  missile  at  the  time  of  impact.  The  process  of  the  filter  is  shown  in  Figures 
3.5,  3.6,  and  3.7. 

3. 5  Frames 

Describing  an  orbit  is  only  possible  with  a  suitable  inertial  reference  frame. 
Coordinate  frames,  as  well  defined  and  detailed  in  many  documents,  differ  from  each 
other  depending  on  their  definition  of  the  origin,  the  fundamental  plane  and  the 
principal  axis.  Only  two  of  these  coordinate  frames,  the  ones  that  are  used  in  this 
document,  will  be  described  in  this  section. 
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Figure  3.6:  Algorithm  for  Least  Squares  Estimation  Filter 

(Part  2  of  3). 


Figure  3.7:  Algorithm  for  Least  Squares  Estimation  Filter 

(Part  3  of  3). 


Body  Frame 


Figure  3.8:  Geocentric-Equatorial  Coordinate  Frame  and 

Satellite  Body  Frame. 

While  describing  the  motion  of  the  objects  orbiting  the  earth,  using  the  geocentric- 
equatorial  coordinate  frame  is  more  convenient  than  other  coordinate  systems.  As 
depicted  in  Figure  3.8,  the  origin  of  the  geocentric-equatorial  coordinate  frame  is  the 
center  of  the  earth  and  its  fundamental  plane  is  the  same  with  the  earth’s  equatorial 
plane.  The  positive  I  axis,  which  is  the  principle  axis,  points  to  the  vernal  equinox 
direction.  The  K  axis  is  perpendicular  to  the  fundamental  plane  and  points  in  the 
direction  of  the  North  Pole.  The  J  axis  completes  the  right  handed  set  of  coordinate 
axes.  Finally,  it  is  important  to  emphasize  that  the  geocentric-equatorial  coordinate 
frame  is  not  rotating  with  the  earth  and  except  for  the  precession  of  the  equinoxes 
the  coordinate  system  is  relatively  fixed  with  respect  to  the  stars. 
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The  body  frame  of  the  satellite  is  used  to  calculate  and  demonstrate  the  position 
of  the  AS  AT  missile  with  respect  to  the  satellite  at  the  time  of  impact.  The  origin  of 
the  body  frame  is  the  center  of  gravity  of  the  satellite.  The  fundamental  plane  of  this 
coordinate  system  is  the  orbital  plate  of  the  satellite.  The  positive  principal  axis  b\  is 
pointing  in  the  direction  of  the  velocity  vector  of  the  satellite.  The  63  axis  is  always 
perpendicular  to  the  orbital  plane.  The  &2  axis  points  roughly  to  the  earth  as  the 
third  element  of  the  coordinate  system,  composed  of  three  perpendicular  unit  vectors. 
As  its  name  implies,  the  body  frame  is  fixed  with  respect  to  the  satellite,  especially 
when  the  two-body  dynamics  assumes  that  the  satellite’s  total  mass  is  concentrated 
at  its  center  of  gravity.  Although  it  is  non-rotating  when  compared  to  the  satellite, 
the  body  frame  rotates  with  respect  to  the  inertial  frame,  stars  and  the  earth.  The 
rotation  matrix  that  converts  a  vector  from  the  initial  frame  to  the  body  frame  can 
be  created  as 


Rbi  =  \ 

w 

{ if 

£ 

■  A 

Z 

■  A 

z 

■  A 

Rbi  = 

Z 

•  A 

A 

■  A 

z 

•  A 

(3.75) 

5 

•  A 

z 

■  A 

Z 1 

•  A 

where  b  =  [61,  62,  A]T  and  i  =  [A,  A,  A]T-  The  vectors  A,  A  and  A  are  the  unit  vectors 
in  the  directions  of  /,  J,  K  respectively. 
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Since  we  have  the  state  vector  of  the  targeted  satellite  in  the  inertial  frame,  the 


principle  axis  of  the  selected  body  frame 
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with  respect  to  inertial  frame  can  be 


found  by  normalizing  the  satellite’s  velocity  vector. 


h 


[VX,  Vy,  Vz 


zisat 


i/BFWW 


(3.76) 


Then,  taking  the  cross  product  of  b\  with  the  position  vector  of  the  satellite  and 

J  i 

normalizing  that  vector,  we  get  [&31  ,  pointing  in  the  direction  of  angular  momentum 


vector, 
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X  [x,y,z]i 
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i 

and  in  order  to  calculate 


cross  product  should  be  used  again 
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(3.78) 


Given 


the  rotation  matrix  from  inertial  to  body  frame  can  be  achieved  easily  in 


two  steps  as 


(3.79) 
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(3.80) 
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The  rotation  matrix  Rbl  was  used  to  express  the  position  of  the  ASAT  on  its  estimated 
trajectory  at  the  time  of  closest  pass. 


- 

“  “ 

X 

X 
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=  Rbi 
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z 

b 

z 

(3.81) 


This  miss  distance  changes  with  the  inputs  of  the  filter.  The  results  will  be  detailed 
in  Chapter  IV. 
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IV.  The  Results  and  the  Analysis 

4-1  Chapter  overview 

In  this  chapter,  the  sensitivity  of  the  filter  was  analyzed  by  comparing  esti¬ 
mated  miss  distances  of  the  ASAT  missile  with  respect  to  the  targeted  satellite.  As 
previously  described,  the  inputs  of  the  filter  that  change  the  estimations  were  the 
sensor  accuracy  and  the  frequency  of  the  observations.  The  filter  was  designed  such 
that  the  accuracy  of  the  sensor  could  be  entered  in  arcseconds  and  the  number  of 
observations  were  entered  as  the  measurement  amount  per  second.  For  each  set  of 
these  two  inputs  the  filter  generated  an  estimated  state  vector  at  epoch  time  and 
this  state  vector  was  propagated  in  time  to  find  the  estimated  state  of  the  ASAT  at 
the  predetermined  impact  time.  This  state  vector  was  defined  and  calculated  with 
respect  to  the  inertial  frame.  Finally  this  vector  was  rotated  and  defined  with  respect 
to  the  body  frame  of  the  satellite.  The  rotation  of  the  state  vector  was  described 


in  the  previous  chapter.  Only  the 


components  of  the  state  vector 


and  the 

were  considered  as  the  cross  track  and  the  radial  components  of  the  miss  distance 
respectively.  For  the  analysis  of  the  results  the  Monte  Carlo  simulation  was  used. 
The  flow  diagram  of  this  analysis  is  shown  in  Figure  4.1. 


Figure  4.1:  Flow  Diagram  of  the  Monte  Carlo  simulations. 
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4-2  Monte  Carlo  Simulations 


Monte  Carlo  simulations  take  the  slightly  different  initial  inputs  for  a  system  and 
iterates  them  in  order  to  find  different  results  based  on  the  accuracy  of  the  iteration 
process  and  dynamics.  After  getting  enough  different  results  some  analysis  like  mean 
and  standard  deviation  of  the  results  can  be  done.  The  Monte  Carlo  simulations  are 
advantageous  because  they  are  easy  to  be  performed  and  only  the  system  dynamics 
are  needed.  [20] 


In  this  research,  Monte  Carlo  simulations  were  used  to  generate  different  esti¬ 
mated  state  vectors  of  the  ASAT  missile  with  respect  to  the  target  satellite  at  the 
impact  time.  At  each  run  of  the  Monte  Carlo  simulation,  the  closest  cross  track  and 
radial  miss  distances  of  the  estimated  ASAT  position  were  calculated.  The  mean 
value  of  the  miss  distances  at  the  plane  perpendicular  to  the  satellite’s  velocity  vector 
was  calculated  as 


%mean 


N 


2=1 


and  standard  deviation  of  the  miss  distances  was  calculated  as 


(4.1) 


&  X 


(Xi  -  x 


)2 

mean  ) 


N 


(4.2) 


where  N  represents  the  amount  of  Monte  Carlo  runs.  In  figures  of  the  results,  an 
ellipsoid,  which  had  an  origin  at  the  mean  of  the  miss  distances,  is  shown.  The  semi¬ 
major  and  the  semi-minor  axis  of  the  ellipsoid  represents  the  standard  deviation  at 
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the  associated  axis  direction.  Also  the  observation  interval  and  the  sensor  accuracy 
is  labeled  in  each  figure.  When  the  input  sets  were  changed  as  more  measurement 
data  per  second  or  as  more  precise  sensor  accuracy  the  miss  distances  became  smaller. 
These  different  input  sets  were  given  to  the  filter  and  after  Monte  Carlo  simulations 
the  outputs  were  calculated  and  shown  in  figures. 

4-3  Comparison  of  Filter  Estimates 

Different  sets  of  inputs  were  created  in  an  observation  amount  range  between  1 
observation  per  second  to  1000  observations  per  second  and  a  sensor  accuracy  range 
between  1  Arcsecod  and  1/1000 Arcsecod.  When  inputs  of  less  than  1  observations  per 
second  was  considered  the  least  squares  filter  was  not  able  to  compute  an  estimate  be¬ 
cause  the  observability  condition  (the  matrix  TtQ~1T  must  be  invertable )  defined  by 
Wiesel  could  not  be  met.  The  sensor  accuracies  more  precise  than  1/1000  Arcsecond 
were  not  considered.  Sensors  with  an  accuracy  of  1  Arcsecond  were  being  used 
at  present  time,  so  anything  better  than  1/1000  Arcsecond  was  considered  beyond 
near-term  capabilities. 

Initially,  the  improvement  of  the  filter  estimates  with  the  sensor  accuracy  based 
on  certain  amount  of  observations  was  described.  Then,  the  comparison  of  the 
estimation’s  improvement  due  to  different  observation  amounts  were  shown  in  this 
section.  For  the  following  different  cases  it  was  seen  that  although  the  filter  had 
a  bias,  the  more  accurate  sensor  or  the  higher  observation  rates  results  in  better 
performance  of  the  filter. 
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4-3.1  1  Observation  per  Second.  For  the  case  in  which  the  the  filter  inputs 

were  composed  of  1  observation  per  second,  it  was  seen  that  the  sensor  accuracy 
should  be  at  least  0.01  Arcsecod.  In  this  case  the  position  estimates  of  the  ASAT 
missile  with  respect  to  the  satellite  were  shown  in  Figure  4.2. 


THE  ASAT  POSITION  WITH  RESPECT  TO  THE  TARGET 


Figure  4.2:  Estimated  position  of  ASAT  at  Impact  Time 

(Data  interval:  lsec;  Accuracy:0.01  Arcsecod). 


When  the  sensor  accuracy  is  increased  to  0.001  Arcseconds  the  miss  distances 
becomes  smaller  as  shown  in  Figure  4.3.  For  these  two  cases  100  Monte  Carlo 
runs  established  and  the  improvement  of  the  filter  estimations  parallel  to  the  sensor 
accuracy  was  observed.  In  reality  the  satellite  would  need  to  stop  taking  observations 
and  make  a  maneuver  to  defeat  the  missile  at  some  time  before  the  impact,  so  the 
estimate  without  this  last  observations  should  be  calculated.  Also,  to  make  the 
scenario  a  little  bit  more  realistic  it  could  be  assumed  that  the  sensor  could  not  take 
measurements  for  some  time  after  the  launch  of  the  missile.  Given  these  facts  the  filter 
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was  tried  to  generate  estimates  with  a  set  of  data  ignoring  some  observations  from 
the  beginning  and  the  last  part  of  the  missile’s  trajectory.  Again  the  observability 
condition  limited  the  filter  and  it  was  seen  that  the  filter  should  have  almost  all  of 
the  observations  in  1  observation  per  second  case. 


THE  ASAT  POSITION  WITH  RESPECT  TO  THE  TARGET 


Figure  4.3:  Estimated  position  of  ASAT  at  Impact  Time 

(Data  interval:  lsec;  Accuracy:!). 001  Arcsecod ). 


The  standard  deviations  of  ASAT  positions,  estimated  at  impact  time  based  on 
one  data  per  second  are  shown  in  Table  4.1. 


Table  4.1:  Standard  Deviations  of  ASAT  Miss  Distances  for 
1  Observation  per  Second. 


Accuracy 

Standard  Deviations  (km) 

(Arcsecond) 

Total 

Cross  Track 

Radial 

0.01 

1.0494 

1.8078 

3.271 

0.001 

0.1203 

0.19987 

0.39567 
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THE  ASAT  POSITION  WITH  RESPECT  TO  THE  TARGET 


Figure  4.4:  Estimated  position  of  ASAT  at  Impact  Time 

(Data  interval:  O.lsec;  Accuracy:!).  1  Arcsecod ). 


4-3.2  10  Observations  per  Second.  When  the  observation  amount  was  in¬ 
creased  it  was  discovered  that  the  filter  could  generate  an  estimate  with  a  senor 
accuracy  of  0.1  Arcsecond.  The  estimated  positions  of  the  ASAT  calculated  with  10 
observations  per  second  and  different  sensor  accuracies  are  shown  in  Figures  4.4,  4.5 
and  4.6. 


It  is  clear  in  the  figures  that  the  estimations  of  the  filter  gets  closer  to  the  refer¬ 
ence  trajectory  and  as  a  result  the  ASAT  missile’s  miss  distance  at  the  impact  time 
gets  smaller  and  smaller  when  the  sensor  becomes  more  accurate.  The  improvement 
of  the  estimates  can  be  seen  together  in  Table  4.2. 
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THE  ASAT  POSITION  WITH  RESPECT  TO  THE  TARGET 


Figure  4.5:  Estimated  position  of  ASAT  at  Impact  Time 

(Data  interval:  O.lsec;  Accuracy:!). 01  Arcsecod ). 


THE  ASAT  POSITION  WITH  RESPECT  TO  THE  TARGET 


Figure  4.6:  Estimated  position  of  ASAT  at  Impact  Time 

(Data  interval:  O.lsec;  Accuracy:0.001  Arcsecod ). 
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Table  4.2:  Standard  Deviations  of  ASAT  Miss  Distances  for 
1 0  Observations  per  Second. 


Accuracy 

Standard  Deviations  (km) 

(Arcsecond) 

Total 

Cross  Track 

Radial 

0.1 

2.7719 

4.3494 

4.316 

0.01 

0.26714 

0.30115 

0.51746 

0.001 

0.03749 

0.060238 

0.037696 

4-3.3  100  Observations  per  Second.  When  the  sensor  begun  to  perform  a 
hundred  obsevations  in  one  second  not  only  the  whole  trajectory  observation  case  but 
also  reduced  data  estimation  case  started  to  establish  better  results.  The  filter  could 


generate  better  estimations  when  data  was  cut  from  the  beginning  and  the  end  of 
the  ASAT  missile’s  trajectory.  The  results  when  the  sensor  could  get  observations 
throughout  the  entire  flight  time  of  the  missile  were  shown  in  Figures  4.7,  4.8,  4.9 
and  the  standard  deviations  of  the  results  for  different  sensor  accuracies  were  shown 


in  Table  4.3. 
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Figure  4.7:  Estimated  position  of  ASAT  at  Impact  Time 

(Data  interval:  O.Olsec;  Accuracy:0. 1  Arcsecod). 
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THE  ASAT  POSITION  WITH  RESPECT  TO  THE  TARGET 


n5 

CL 


Figure  4.8:  Estimated  position  of  ASAT  at  Impact  Time 

(Data  interval:  O.Olsec;  Accuracy:!). 01  Arcsecod ). 


-5.  THE  ASAT  POSITION  WITH  RESPECT  TO  THE  TARGET 


Figure  4.9:  Estimated  position  of  ASAT  at  Impact  Time 

(Data  interval:  O.Olsec;  Accuracy:0.001  Arcsecod ). 
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Table  4.3:  Standard  Deviations  of  ASAT  Miss  Distances  for 
1 00  Observations  per  Second. 


Accuracy 

Standard  Deviations  (km) 

(Arcsecond) 

Total 

Cross  Track 

Radial 

0.1 

0.68924 

1.2357 

0.32363 

0.01 

0.080419 

0.18342 

0.027629 

0.001 

0.015154 

0.015154 

0.003482 

4-3-4  1000  Observations  per  Second.  When  the  case  in  which  the  sensor 

performs  1000  observations  per  second  was  considered  the  total  amount  of  measure¬ 
ments  for  this  particular  scenario  reached  up  to  500,  000.  Since  the  length  of  the 
matrixes  was  too  large  the  numerical  integration  became  difficult  to  handle  with  or¬ 
dinary  home-use  computers  with  limited  processors  and  memory.  The  problem  was 
still  solvable  but  took  a  lot  of  time  (i.e.  approximately  two  days)  to  get  the  results 
of  the  filter’s  one  or  two  estimates.  Because  of  the  limited  time  of  this  research  only 
two  or  three  filter  estimates  based  on  input  sets  of  1000  observations  per  second  and 
different  sensor  accuracies  were  generated.  The  results  were  shown  in  Figures  4.10, 
4.11,  4.12.  Even  though  only  two  or  three  computations  could  be  established  the 
results  were  shown  in  the  figures  and  the  standard  deviations  for  each  axis  had  been 
computed  and  summarized  in  Table  4.4  in  order  to  compare  with  the  other  cases. 
The  standard  deviations  could  not  be  accurate  enough  due  to  the  limited  runs  of  the 
Monte  Carlo  simulation.  However,  the  order  of  the  miss  distances  of  ASAT  missile 
at  impact  time  can  give  a  sense  and  make  the  comparisons  more  meaningful.  The 
mean  values  of  the  miss  distances  are  shown  in  Table  4.5. 
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THE  ASAT  POSITION  WITH  RESPECT  TO  THE  TARGET 


Figure  4.10:  Estimated  position  of  ASAT  at 

(Data  interval:  O.OOlsec;  Accuracy:!  Arcsecod ). 


Impact  Time 


Cross  Track  Misdistance  (Km.)  [b2 - >(Earth)] 


Figure  4.11: 
(Data  interval: 


Estimated  position  of  ASAT  at  Impact  Time 
O.OOlsec;  Accuracy:0.1  Arcsecod ). 
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THE  ASAT  POSITION  WITH  RESPECT  TO  THE  TARGET 


Figure  4.12:  Estimated  position  of  ASAT  at  Impact  Time 

(Data  interval:  O.OOlsec;  Accuracy:!). 01  Arcsecod ). 


Table  4.4:  Standard  Deviations  of  ASAT  Miss  Distances  for 
1000  Observations  per  Second. 


Accuracy 

Standard  Deviations  (km) 

(Arcsecond) 

Total 

Cross  Track 

Radial 

1 

0.69424 

3.2288 

0.012505 

0.1 

0.12873 

0.50187 

0.00809 

0.01 

0.0326 

0.0326 

0.00045 

Table  4.5:  Mean  Values  of  ASAT  Miss  Distances  for  1000 

Observations  per  Second. 


Accuracy 

Mean  Values  (km) 

(Arcsecond) 

Total 

Cross  Track 

Radial 

1 

3.23 

0.69 

0.13 

0.1 

0.51 

0.16 

0.064 

0.01 

0.061 

0.0612 

0.066 

4-3.5  The  Filter  Estimates  in  More  Realistic  Simulations  .  As  mentioned 
previously,  the  satellite  should  stop  taking  measurements  at  some  time  before  the 
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predetermined  impact  time  in  order  to  maneuver  and  defeat  the  ASAT  missile.  Also, 
in  a  more  realistic  simulation  the  satellite  would  only  be  able  to  catch  and  start 
tracking  the  ASAT  missile  at  some  time  after  the  launch.  The  filter  was  tested  for 
different  cases  in  which  the  observation  started  at  a  time  after  the  launch  and  the 
measurements  stopped  at  a  time  before  the  impact.  In  most  trials  the  observability 
condition  could  not  be  met  and  the  filter  could  not  generate  an  estimated  ASAT 
trajectory.  It  has  been  concluded  that  as  an  input  to  the  filter,  the  observation 
amount  should  be  at  least  10  observations  per  second  in  order  to  accomplish  an 
estimate  without  observing  the  entire  trajectory  from  the  launch  to  the  impact.  The 
different  input  sets  that  resulted  with  an  estimate  were  pictured  in  the  following 
figures.  For  the  first  condition  the  filter  could  calculate  an  estimated  trajectory 
by  tracking  the  missile  from  60  seconds  after  the  launch  until  10  seconds  before  the 
impact  time  with  an  observation  rate  of  100  observations  per  second. 

In  this  case,  if  we  assume  that  the  satellite  could  start  the  maneuver  just  after  it 
stopped  taking  the  observations,  again  assuming  that  it  would  maneuver  in  the  radial 
direction  where  the  standard  deviation  of  the  estimated  miss  distances  is  larger  than 
the  cross  track  direction  and  with  an  assumption  that  it  would  maneuver  just  enough 
to  move  out  of  the  standard  deviation  amount  of  range, 
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THE  ASAT  POSITION  WITH  RESPECT  TO  THE  TARGET 


Figure  4.13:  (Data  interval:  O.Olsec; 

Accuracy:!).  1  Ar cs e cod: from  60  sec  after  launch  till  10  sec 
to  impact). 


d  =  V0 1+  ^ gt 2  (4.3) 

5000m  =  ^(10)2 
g  «  100  m/ sec2 


the  acceleration  on  the  satellite  will  be  approximately  100  m/sec2.  In  order  to  generate 


this  acceleration  the  thrusters  on  a  satellite,  which  is  in  the  size  of  FY-1C,  should 


create  the  force, 
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F  =  ma  (4.4) 

«  750  kg  x  100  ml  sec? 

F  «  75000  IV 

Although  this  calculated  thrust  seems  to  be  a  little  bit  high,  there  are  available  engines 
that  can  create  it.  [21]  It  can  be  said  that  if  the  accuracy  of  the  sensor  would  be 
increased,  more  precise  estimations  can  be  generated  by  the  filter,  but  since  this  result 
is  reasonable  and  applicable,  the  more  accurate  cases  wasn’t  calculated. 

For  the  next  case  if  the  observation  amount  was  reduced  to  ten  measurements 
per  second  then  the  sensor  accuracy  should  be  increased  in  order  to  have  the  filter 
generate  an  estimate.  The  calculated  estimates  were  pictured  in  Figures  4.14,  4.15, 
4.16.  Using  the  same  maneuver  assumptions  made  for  the  previous  case  and  the 
Equations  4.3  and4.4,  if  the  sensor  performed  10  observations  per  second  and  started 
to  track  the  ASAT  missile  60 seconds  after  the  launch;  the  approximate  amount  of 
maneuver  accelerations  of  the  satellite  and  the  required  forces  would  be  like  the  cal¬ 
culated  results  shown  in  Table  4.6.  The  required  thrusts  that  for  the  associated 
accelerations  were  calculated  based  on  a  750  kg  satellite  representing  FY-1C. 


79 


THE  ASAT  POSITION  WITH  RESPECT  TO  THE  TARGET 


Figure  4.14:  (Data  interval:  O.lsec; 

Accuracy:!). 01  Arcsecod] from  60  sec  after  launch  till  10 
sec  to  impact). 


THE  ASAT  POSITION  WITH  RESPECT  TO  THE  TARGET 


Figure  4.15:  (Data  interval:  O.lsec; 

Accuracy:0.001  Arcsecod] from  60  sec  after  launch  till  10 
sec  to  impact). 
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THE  ASAT  POSITION  WITH  RESPECT  TO  THE  TARGET 


Figure  4.16:  (Data  interval:  O.lsec; 

Accuracy:!). 001  Arcsecod/.irom  60  sec  after  launch  till  20 
sec  to  impact). 


Table  4.6:  Required  Accelerations  and  Forces  of  Satellite 

Maneuvers  for  10  observations  per  Second. 


Accuracy 

First  Data 

Last  Data 

Aceleration 

Force 

Arcsecond 

sec  after  Launch 

sec  before  Impact 

m / sec2 

N 

0.01 

60 

10 

~  40 

«  30000 

0.001 

60 

10 

«  4.4 

«  3300 

0.001 

60 

20 

~  1.4 

~  1050 

4-3.6  Comparison  Based  on  Data  Interval  or  Sensor  Accuracy  Only.  In 
order  to  compare  the  estimates  of  the  least  squares  filter,  the  results  can  be  compared 
while  keeping  the  accuracy  constant  and  changing  the  observation  amount  per  second. 
The  improvement  of  the  standard  deviations  can  provide  an  idea  about  the  sensor 
selection.  Since  observability  condition  could  not  be  met  at  each  set  of  the  inputs, 
only  the  results  of  the  achievable  input  sets  were  shown.  For  the  first  condition 
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when  the  accuracy  was  0.1  Arcsecond  the  standard  deviation  of  the  miss  distances  got 


smaller  while  the  observation  amount  increased,  as  seen  in  Figure  4.17. 


Improvement  of  Standard  Deviations  (0.1  Arcsecond  Accuracy) 


Figure  4.17:  Improvement  of  Standard  Deviation  (0.1  Arcsec¬ 
ond  Accuracy). 

When  the  accuracy  was  increased  to  0.01  Arcsecond  the  standard  deviation  improve¬ 
ment  got  better  as  it  could  be  seen  in  Figure  4.18.  Figure  4.19  shows  the  case  in 
which  the  accuracy  had  been  improved  to  0.001  Arcsecods. 

Finally,  with  the  results  that  could  be  calculated,  one  more  comparison  could  be 
done  by  keeping  the  observation  amount  per  second  constant  and  changing  the  accu¬ 
racy.  These  comparisons  are  shown  in  the  Figures  4.20  and  4.21.  It  was  evident  in  the 
figures  that  even  though  the  improvement  in  the  sensor  accuracy  effected  the  filter’s 
estimations  considerably,  the  effect  of  accuracy  growth  reduced  after  0.01  Arcsecond , 
giving  a  sense  that  increasing  accuracy  further  from  this  point  would  not  help  much. 
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Improvement  of  Standard  Deviations  (0.01  Arcsecond  Accuracy) 


Figure  4.18:  Improvement  of  Standard  Deviation  (0.01  Arc- 
second  Accuracy). 


Improvement  of  Standard  Deviations  (0.001  Arcsecond  Accuracy) 


Figure  4.19:  Improvement  of  Standard  Deviation  (0.001  Arc- 
second  Accuracy). 
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Improvement  of  Standard  Deviations  (10  observations  per  Second) 


Figure  4.20:  Improvement  of  Standard  Deviation  (10  Obser¬ 
vations  per  Second). 


Improvement  of  Standard  Deviations  (10  observations  per  Second) 


Figure  4.21:  Improvement  of  Standard  Deviation  (100  Obser¬ 
vations  per  Second). 
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V.  Conclusions 


5. 1  Summary 

This  thesis  demonstrates  the  ability  to  estimate  the  trajectory  of  a  threat  ASAT 
missile  based  on  the  angle-only  observations  conducted  by  a  space-based  onboard 
optical  sensor.  The  importance  of  having  a  self  defense  capability  for  a  satellite 
was  emphasized  with  a  brief  description  of  the  background  of  ASAT  systems,  space 
surveillance  systems  and  the  previous  studies  about  this  topic.  This  research  aims 
to  equip  a  space  asset  with  an  ability  to  first  detect  and  track  an  ASAT  missile  and 
then  estimate  its  future  position  with  respect  to  its  own  position.  Based  on  this 
capability  the  space  asset  would  have  a  better  situational  awareness  and  probably 
have  time  for  a  last  ditch  maneuver  to  defeat  the  ASAT  missile  and  survive.  A  least 
squares  estimation  filter  was  created  and  tested  on  a  particular  simulated  scenario  for 
different  sensor  specifications  in  order  to  determine  the  required  sensor’s  capabilities. 

5.2  Conclusions 

The  estimates  of  the  filter  were  analyzed  in  the  previous  chapter.  Since  the 
observation  inputs  of  the  designed  filter  were  provided  by  only  one  optical  sensor  and 
the  sensor  was  assumed  to  perform  angle-only  measurements,  it  has  been  observed 
that  the  filter  was  highly  limited  due  to  the  observability  condition,  as  described  in 
Chapter  III .  Simulation  studies  showed  that  the  sensor  accuracy  should  be  at  least 
1  Arcsecond  to  be  able  to  get  an  estimate  from  the  filter  which  gets  a  maximum 
observation  amount  of  1000  observations  per  second.  From  another  point  of  view  the 
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sensor  should  perform  at  least  one  observations  per  second  for  the  proper  processing  of 
this  estimation  filter  using  a  sensor  accuracy  of  maximum  1/1000  Arcsecond.  These 
conclusions  were  for  the  case  in  which  the  sensor  could  perform  observations  from  the 
very  beginning  until  the  impact  point  of  the  ASAT’s  trajectory. 

To  continue  with  the  whole  observation  of  the  ASAT’s  trajectory  case,  the 
analysis  of  the  simulations  shows  that  if  the  sensor  can  perform  10  observations  per 
second  with  an  accuracy  of  0.1  Arcsecond  the  standard  deviations  of  the  estimated 
miss  distances  in  the  radial  axis  (worst  case)  will  be  approximately  4.5  km.  The 
satellite  can  maneuver  out  of  this  danger  zone  with  an  acceptable  acceleration  rate. 
Also  it  was  observed  that  increasing  the  observation  rate  or  the  accuracy  of  the 
sensor  beyond  these  values  did  not  improve  the  results  considerably.  But  with  the 
consideration  that  the  observations  can  be  increased  up  to  100  hundred  per  second 
easily,  this  range  can  be  reduced  down  to  approximately  1.3 km  with  using  a  sensor 
with  the  same  accuracy  of  0.1  Arcsecond. 

Although  the  filter  needs  continuous  observation  of  the  ASAT  in  most  cases, 
some  more  realistic  cases  could  be  generated  to  evaluate  the  capability  of  the  filter 
and  compare  the  possible  sensor  options.  The  main  objective  in  these  simulations  was 
to  leave  the  satellite  some  time  before  the  impact  to  maneuver  and  also  to  consider 
the  fact  that  the  sensor  can  start  the  observations  at  some  time  after  the  launch.  The 
simulations  showed  that  the  sensor  accuracy  should  be  at  least  0.1  Arcsecond  in  order 
to  be  able  to  reduce  observation  data  and  have  the  filter  generate  an  estimate  within 
the  observation  amount  and  the  sensor  accuracy  range  that  have  been  tested.  Using 


a  0.1  Arcsecond  accurate  sensor,  the  analysis  showed  that  if  it  was  assumed  that  the 
sensor  could  start  performing  observations  1  minute  after  the  launch  and  stopped  the 
observations  10  seconds  before  the  impact  time,  with  100  observations  per  second  the 
standard  deviation  of  the  estimated  miss  distances  became  ~  5  km.  Assuming  that 
the  satellite  could  start  the  maneuver  immediately,  the  defeat  maneuver  could  be  done 
with  an  acceleration  rate  of  ~  lOg  's.  If  the  observation  frequency  was  reduced  to 
10  observations  per  second  and  sensor  accuracy  was  increased  to  be  0.01  Arcsecond 
then  the  required  acceleration  rate  became  ~  4g  's. 

Finally,  it  was  concluded  that  a  space-based  onboard  sensor  with  an  accuracy  of 
0.1  Arcsecond  and  the  capability  to  perform  100  observations  per  second  can  establish 
the  proposed  operation  for  this  specific  scenario  with  reasonable  results. 

5. 3  Future  Work 

During  the  simulations  it  has  been  discovered  that  the  observability  condition 
limits  the  filter  estimates.  If  another  observation  can  be  added  to  the  calculations  the 
observability  range  can  be  enlarged.  The  research  results  can  be  improved  by  using 
the  same  filter  and  additional  observation  measurements.  Additional  observations 
can  be  performed  by  a  sensor  mounted  on  a  geostationary  satellite,  mounted  on  a 
satellite  in  the  same  formation  (if  there  is  one)  or  another  sensor  on  the  same  satellite 
with  a  lateral  distance.  Also,  the  effect  of  using  an  active  optical  sensor  and  a  radar 
sensor  on  the  estimations  should  be  considered  and  calculated. 
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The  perturbations  on  the  two-body  motion  were  not  calculated  in  this  research, 
so  in  a  future  work  the  effect  of  the  perturbations  and  atmospheric  effects  should  be 
included  into  the  calculations.  Also  the  observations  obstructed  by  the  curvature  of 
the  earth  should  be  considered. 

The  simulations  in  this  research  were  conducted  using  ordinary  home-use  com¬ 
puters.  Based  on  this  fact  only  a  limited  amount  of  simulations  could  be  performed 
in  a  limited  time.  More  sets  of  inputs  can  be  generated,  more  simulations  can  be  per¬ 
formed  and  the  requirements  of  the  optimum  sensor  can  be  determined  in  a  narrower 
range  if  computers  with  faster  processors  and  larger  memories  are  used. 

The  non-linear  least  squares  filter  produced  in  this  research  is  a  batch  filter  and 
it  has  to  wait  until  all  of  the  data  is  available  to  make  the  estimate.  In  real  life,  the 
satellite  would  want  to  continuously  update  the  estimate  of  the  ASAT’s  trajectory. 
That  would  lead  to  a  sequential  type  of  filter  like  a  Kalman  Filter.  As  a  start  in  this 
research  topic,  using  a  more  stable  filter  like  least  squares  to  figure  out  the  quality 
and  the  quantity  requirements  of  the  sensor  was  considered  more  appropriate.  Given 
the  recommendations  for  the  quality  and  the  quantity  of  the  data,  a  possible  future 
work  will  have  a  good  starting  point  for  working  with  the  Kalman  filter. 

In  real-time  processing,  it  is  important  for  the  computer  to  “out  run”  real  life. 
Also,  in  this  scenario  the  faster  the  data  comes  in,  the  less  time  there  is  for  the 
on-board  computer  to  calculate  the  ASAT’s  trajectory.  This  is  another  reason  to 
apply  a  sequential  filter  into  these  calculations.  While  working  on  the  sequential 


filtering,  a  future  research  should  also  examine  whether  a  typical  computer  can  “out 
run”  realtime  or  not.  [10] 

This  research  examined  only  the  case  in  which  the  ASAT  was  assumed  to  be 
an  unguided  ballistic  missile.  Precautions  against  other  kind  of  ASAT  systems  like 
LASER  ASATs  or  guided  KE  ASATs  should  be  considered  in  the  future  studies. 


Appendix  A.  MATLAB  Code 


A .  1  Main  Code 


clc, clear, elf; 
load  data 
du=6378. 135; 

tu=sqrt (6378 . 135~3/3 . 98601e5) ; 
mu=l ; 


7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7 

/o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /  o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /  o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 

datapersec=100 ;  Renter  observation  amount  per  second 
first=60  Renter  fist  second  of  data 

last=10  Renter  last  data  second  (how  many  sec  before  impact) 

accuracy= .  001  °/0arcseconds 

compute_z(datapersec) ; 


0/  0/  0/  0/  0/  0/ 


0  / 0  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  / 0  /O  /O  /O  /O  /O  /O  /O  /O  / 0  /O  /O  /O  / 0  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 


f irstdata=f irst*datapersec 

if  datapersec==l 
load  zQdata_lsec 

z_i=z_i_lsec_p001 (f irstdata: length(z_lsec) , : ) ; 
Q=Q_lsec_p001 (f irstdata: length(z_lsec) , ; 
lastdata=length(z_lsec)-f irstdata- last *datapersec 
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end 


if  datapersec==10 
load  zQdata_l 

z_i=z_i_l_p001 (f irstdata: length(asatl) , ; 
Q=Q_l_p001 (f irstdata: length(asatl) , : ) ; 
lastdata=length(asatl) -f irstdata-last*datapersec 
end 

if  datapersec==25 
load  zQdata_p04sec 

z_i=z_i_p04sec_pl (f irstdata: length(z_p04sec) , : ) ; 
Q=Q_p04sec_pl (f irstdata: length(z_p04sec) , : ) ; 
lastdata=length(z_p04sec) -f irstdata-last*datapersec 
end 

if  datapersec==100 
load  zQdata_01 

z_i=z_i_01_p001 (f irstdata: length (asat02) ,  : )  ; 
Q=Q_01_p001 (f irstdata: length(asat02) , : ) ; 
lastdata=length(asat02)-f irstdata-last*datapersec 
end 
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for  clmn=l:10 


[state_last , target] =statelast (f irstdata, lastdata,z_i , Q , clmn,datapersec) ; 

clmn=clmn 

hold  on 

[pos ,pos_cross ,pos_radial] =compute_r (state_last , clmn, target) ; 
hold  on 

p_all (clran)=pos ; 
cross (clmn)=pos_cross ; 
radial (clmn) =pos_radial ; 

end 


[mean, stdev]  =  stat(p_all); 

[mean_cross , stdev_cross]  =  stat (cross); 
[mean_radial , stdev_radial]  =  stat (radial) ; 

hold  on 
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ellipsedraw(stdev_pl_pl_cross*du, stdev_pl_pl_radial*du, . . . 
mean_pl_pl_cross*du,mean_pl_pl_radiai*du,0, ’ r’ ) ; 

hold  on 

xlabel(’ Cross  Track  Misdistance  (Km.)  [b2 - > (Earth) ]’) 

ylabel ( ’Radial  Misdistance  (Km.)  [b3 - >(normal  to  orbital  plane)]’) 

grid  on 

title ( ’THE  ASAT  POSITION  WITH  RESPECT  TO  THE  TARGET’ , ’FontSize’ , 10) 
str(l)  =  {[’Data  interval :’ ,num2str (1/datapersec) ,  ’ sec ’]} ; 
str(2)={ [’SensorAccuracy: ’ ,num2str (accuracy) , ’Arcsec’] }; 
str(3)={ [’FirstData: ’ ,num2str(f irst) , ’sec  After  Launch’]}; 
str(4)={ [’LastData: ’ ,num2str (last) , ’  sec  to  impact’]}; 
str(5)={ [’Std.Dev. : ’ ,num2str (stdev*du) , ’  Km’] }; 
str(6)={ [’Std.Dev.CrossTrack: ’ ,num2str (stdev_cross*du) , ’Km’] }; 
str (7)={ [’Std.Dev . Radial : ’ ,num2str (stdev_radial*du) , ’Km’] } ; 
str(8)  =  {’Ellipse  :  Std.Dev. -center  at  mean’}; 
h  =  axes ( ’Position’ ,  [0  0  1  1] , ’Visible ’,’ of f ’) ; 
set (gcf , ’ Current Axes ’ , h) 
text ( . 65 , . 79 , str, ’FontSize ’ ,7) 
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A. 2  MATLAB  Functions  Used  in  the  Main  Code 


A. 2.1  The  Function  to  Create  the  Observations. 

function  [z_all] =compute_z (datapersec) 

du=6378. 135; 

tu=sqrt (6378 . 135~3/3 . 98601e5) ; 
mu=l ; 
load  data 

/T/„  1/10  SECOND 

if  datapersec==10 

asat= [asatl( : ,4)/du  asatl ( : , 5) /du  asatl ( : , 6)/du] ; 
t= [0 : . 1/tu: ( ( (asatl (length (asat) , 1) -asatl (1 , 1) ) *24*3600)+ . 1) /tu] ; 
target= [target 1 ( : , 4)/du  target 1 ( : , 5) /du  target 1 ( : , 6)/du  ]; 
z_all=zeros (length(asat) , 100) ; 

Q_all=zeros (length(asat) , 100) ; 
for  n= [1  .1  .01  .001] ; 

sigma_sens=n/3600*pi/ 180 ; 
c=random( ' norm ’ , 0 , sigma_sens , length (asat) , 100) ; 
for  k=l : length(asat) 
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z(k, l)=abs ( (asat (k, 1 : 3) ) ^target (k, 1 : 3) 5 /(norm(asat (k, 1 : 3) ) * . 


norm(target (k, 1 : 3) ) ) ) ; 


for  m=l : 100 

z_i (k,m)=z(k, l)+c (k,m)*sqrt (1-z (k, 1) "2) ; 
if  z_i(k,m)>l 

z_i (k,m)=z(k, 1) ; 

end 

Q (k,m)=(sigma_sens~2) * (1- (z_i (k,m) ) ~2) ; 

end 

end 

z_all=[z_all  z_i] ; 

Q_all= [Q_all  Q] ; 
a=length (z_all ( 1 , : ) ) 
end 

z_l=z; 

z_i_l_l=z_all ( : ,101:200)  ; 
z_i_l_pl=z_all ( : ,201 :300)  ; 
z_i_l_p01=z_all ( : , 301 : 400)  ; 
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z_i_l_p001=z_all ( : ,401:500) ; 

Q_l_l=Q_all(: ,101:200) ; 

Q_l_pl=Q_all(: ,201:300) ; 

Q_l_p01=Q_all(: ,301:400) ; 

Q_l_p001=Q_all(: ,401:500) ; 

clear  z_all  Q_all  a  target  t  asat  asatl  asat02... 

asat002  target  1  target02  target002  c  k  mu  du  tu  sigma_sens  m  n  z 
save  zQdata_l 
end 


111  1/100  SECOND 

if  datapersec==100 
k=l ; 

asat= [asat02( : ,4)/du  asat02( : ,5)/du  asat02( : ,6)/du] ; 
target= [target02 ( : ,4)/du  target02( : ,5)/du  target02( : ,6)/du] ; 
z_all=zeros (length(asat) , 10) ; 

Q_all=zeros (length(asat) , 10) ; 
for  n=  [1  10  100  .1  .01  .001] ; 
sigma_sens=n/3600*pi/ 180 ; 

c=random( 'norm’ , 0 , sigma_sens , length (asat) , 10) ; 
for  k=l : length (asat) 
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z(k, l)=(asat (k, 1 : 3) ) *target (k, 1 : 3) V (norm(asat (k, 1 : 3) ) * . 


norm (target (k , 1 : 3) ) ) ; 


for  m=l : 10 

z_i (k,m)=z(k, l)+c (k,m) *sqrt (l-z(k, 1) ~2) ; 
if  z_i(k,m)>l 

z_i(k,m)=z(k, 1) ; 

end 

Q (k,m)=(sigma_sens~2) * (1- (z_i (k,m) ) ~2) ; 

end 

end 

z_all=[z_all  z_i] ; 

Q_all=[q_all  q] ; 
a=length (z_all ( 1 , : ) ) 

end 

z_01=z ; 

z_i_01_l=z_all ( : ,11:20)  ; 
z_i_01_10=z_all ( : ,21:30) ; 
z_i_01_100=z_all ( : ,31:40) ; 
z_i_01_pl=z_all ( : ,41 :50) ; 
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z_i_01_p01=z_all ( : ,51 : 60) ; 


z_i_01_p001=z_all ( : , 61 : 70) ; 

Q_01_l=Q_all ( : ,11:20) ; 

Q_01_10=Q_all(: ,21:30) ; 

Q_01_100=Q_all(: ,31:40) ; 

Q_01_pl=Q_all(: ,41:50) ; 

Q_01_p01=Q_all(: ,51:60)  ; 
q_01_p001=q_all ( : , 61 : 70) ; 

clear  z_all  q_all  a  target  t  asat  asatl  asat02  asat002... 

targetl  target02  target002  c  k  mu  du  tu  sigma_sens  m  n  z 
save  zqdata_01 

end 


du=6378. 135; 

tu=sqrt (6378 . 135~3/3 . 98601e5) ; 
mu=l ; 


load  data 


°/„°/o7o  1  SECOND 

if  datapersec==l 
for  j=0:524; 

asat (j+1 , 1 : 6)= [asatl (10*j+l ,4) /du  asatl(10*j+l,5)/du. . . 

asatl(10*j+l,6)/du  asatl (10*j+l , 7) *tu/du  asatl (10*j+l , 8) *tu/du. 
asatl(10*j+l,9)*tu/du] ; 

target (j+1 , 1 : 6)= [targetl (10*j+l ,4) /du  targetl(10*j+l,5)/du. . . 
targetl (10*j+l , 6) /du  targetl(10*j+l,7)*tu/du. . . 
targetl (10*j+l , 8) *tu/du  targetl (10* j+1 , 9) *tu/du] ; 

end 

t=  [0 : 1/tu: j/tu] ; 
z_all=zeros (length(asat) , 100) ; 

Q_all=zeros (length(asat) , 100) ; 
for  n= [1  .1  .01  .001] ; 


sigma_sens=n/3600*pi/ 180 ; 

c=random( 'norm’ , 0 , sigma_sens , length (asat) , 100) ; 
for  k=l : length (asat) 

z(k, l)=abs ( (asat (k, 1 : 3) ) ^target (k, 1 : 3) 5 /(norm (asat (k, 1 : 3) ) * . . . 


norm(target (k, 1 : 3) ) ) ) ; 


for  m=l : 100 


z_i (k,m)=z(k, l)+c(k,m) *sqrt (l-z(k, 1)  ~2) 
if  z_i(k,m)>l 

z_i(k,m)=z(k, 1) ; 

end 

Q(k,m)=(sigma_sens~2) * (l-(z_i (k,m) ) ~2) ; 

end 


end 

z_all=[z_all  z_i] ; 

Q_all= [Q_all  Q] ; 
a=length (z_all ( 1 , : ) ) 
end 

z_lsec=z ; 

z_i_lsec_l=z_all ( : , 101 : 200) ; 
z_i_lsec_pl=z_all ( : , 201 : 300) ; 
z_i_lsec_p01=z_all ( : ,301 :400) ; 
z_i_lsec_p001=z_all( : ,401 : 500) ; 
clear  z_all 
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Q_lsec_l=Q_all ( : ,101:200) ; 


Q_lsec_pl=Q_all(: ,201:300) ; 

Q_lsec_p01=Q_all(: ,301:400) ; 

Q_lsec_p001=Q_all ( : ,401:500) ; 

clear  z_all  Q_all  a  target  t  asat  asatl  asat02  asat002... 

targetl  target02  target002  c  k  mu  du  tu  sigma_sens  m  n  z 
save  zQdata_lsec 
end 


A. 2.2  The  Function  to  Compute  the  Estimated  State  at  the  Impact  Time. 

function  [state_last,target]=statelast(firstdata,lastdata,z_i, . . . 

Q , clmn , datapersec) 

delta_x=le20* [1 ; 1 ; 1 ; 1 ; 1 ; 1] ; 

P=eye (6) ; 

%using  canonical  units 
du=6378. 135; 

tu=sqrt (6378 . 135~3/3 . 98601e5) ; 
mu=l ; 
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load  data 


HI  1  SECOND 
if  datapersec==l 
for  j=0:524; 

asat (j+1 , 1 : 6)= [asatl (10*j+l ,4) /du  asatl(10*j+l,5)/du. . . 
asatl (10*j+l , 6)/du  ... 

asatl(10*j+l,7)*tu/du  asatl(10*j+l,8)*tu/du. . . 
asat 1 ( 10* j +1 , 9) *tu/ du] ; 

target (j+1 , 1 : 6)  =  [target  1 (10* j+1 ,4) /du  targetl(10*j  +  l,5)/du. . . 
targetl (10*j+l , 6)/du  targetl(10*j+l,7)*tu/du. . . 
targetl (10*j+l , 8) *tu/du  targetl (10*j+l , 9) *tu/du] ; 
end 

target= [target (firstdata: length(asat) , 1) . . . 
target (firstdata: length (asat) ,2) . . . 
target (firstdata: length (asat) ,3)  . .  . 
target (firstdata: length (asat) ,4) . . . 
target (firstdata: length (asat) ,5) . . . 
target (firstdata: length (asat) ,6)] ; 
asat= [asat (firstdata: length (asat) , 1) . . . 
asat (firstdata: length (asat) ,2) . . . 
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asat(f irstdata: length (asat) ,3) 


asat(f irstdata: length(asat) ,4) . . . 
asat (f irstdata: length(asat) ,5) . . . 

asat (f irstdata: length (asat) ,6)] ; 
t=  [0 : 1/tu: (length(asat)-l)/tu] ; 

state_zero= [asat (1 , 1)  asat(l,2)  asat(l,3)  asat (1,4)  asat(l,5)  asat(l,6)]5; 

end 


T/l  1/10  SECOND 
if  datapersec==10 

asat= [asatl (f irstdata: length(asatl) ,4)/du  . . . 
asatl (f irstdata: length(asatl) ,5)/du  ... 
asatl (f irstdata: length (asatl) ,6)/du  ... 
asatl (f irstdata: length (asatl) ,7)*tu/du. . . 
asatl (f irstdata: length (asatl) , 8)*tu/du. . . 
asatl (f irstdata: length (asatl) , 9) *tu/du] ; 
t=  [0 : . 1/tu: ( ( (asatl (length (asatl) , 1)- . . . 
asatl (f irstdata, 1) ) *24*3600)+ . l)/tu] ; 
target= [targetl (f irstdata: length(asatl) ,4)/du. . . 
target 1 (f irstdata: length (asatl) , 5)/du. . . 
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target 1 (firstdat a: length (asatl) , 6)/du  ... 
target 1 (firstdat a: length (asatl) ,7) *tu/du. . . 
targetl (f irstdata: length(asatl) ,8) *tu/du. . . 

targetl (f irstdata: length(asatl) , 9) *tu/du] ; 
state_zero= [asat (1 , 1)  asat(l,2)  asat(l,3)  asat(l,4)  asat(l,5)  asat(l,6)] ’ 

end 


111  1/100  SECOND 
if  datapersec==100 

asat= [asat02(f irstdata: length (asat02) ,4) /du. . . 
asat02(f irstdata: length (asat02) ,5)/du. . . 
asat02(f irstdata: length (asat02) ,6)/du  ... 
asat02(f irstdata: length (asat02) , 7)*tu/du. . . 
asat02(f irstdata: length (asat02) ,8)*tu/du. . . 
asat02(f irstdata: length (asat02) , 9)*tu/du] ; 
t=  [0 : . 01/tu: ( ( (asat02(length(asat02) , 1)  - .  .  . 
asat02(f irstdata, l))*24*3600)+.01)/tu] ; 
target= [target02 (f irstdata : length(asat02) ,4)/du. . . 


target02 (f irstdata: length(asat02) ,5)/du. . . 


target02 (f irstdata: length(asat02) ,6)/du  . . . 
target02 (f irstdata: length(asat02) ,7)*tu/du. . . 
target02(f irstdata : length(asat02) ,8)*tu/du. . . 
target02 (f irstdata: length(asat02) ,9)*tu/du] ; 

state_zero= [asat (1 , 1)  asat(l,2)  asat(l,3)  asat(l,4)  asat(l,5)  asat(l,6)]’; 
end 


0/  0/  0/  0/  0/  0/ 


0  /0  /O  /O  / 0  / 0  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  / 0  /O  /O  /O  /O  /O  /O  /O  /O  / 0  /O  /O  /O  /O  /O  /O  /O  /O  / 0  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  /O  / 0  /O  /O  /O  /O  /O  /O  /O  /O  / 0  /o  /o  /o  /o  /o  /o  /o  /o  /o  /o 


state_new=state_zero ; 


while  length (nonzeros (delta_x< . l*sqrt (diag (P) ) ) ) <6 


state_zero=state_new 

phi_zero=[l  zeros(l,6)  1  zeros(l,6)  1  zeros(l,6)  1  zeros(l,6)  1... 
zeros(l,6)  1] ; ; 

y_zero= [phi_zero ; state_zero]  ’  ; 

[t ,y] =ode45(@y_eom_l ,t ,y_zero)  ; 


105 


P_inv=0; 


residual=0 ; 
rsdl_sum=0 ; 
delta_x=0 ; 

for  i=(l : length(t) ) 

state (i, l:6)=y(i,37:42) ; 

phi=  [y(i , 1)  y (i , 2)  y(i,3)  y(i,4)  y(i,5)  y (i , 6) ; . . . 
y (i , 7)  y(i,8)  y(i,9)  y(i,10)  y(i,ll)  y (i , 12) ; . . . 
y (i , 13)  y (i , 14)  y(i,15)  y(i,16)  y(i,17)  y (i , 18) ; . . . 
y (i , 19)  y (i , 20)  y(i,21)  y(i,22)  y(i,23)  y(i,24);... 
y(i ,25)  y (i , 26)  y(i,27)  y(i,28)  y(i,29)  y(i,30);... 
y (i , 31)  y (i , 32)  y(i,33)  y(i,34)  y(i,35)  y (i , 36) ; ]  ; 

denum=sqrt  ( (stated ,  1)  “2+state(i  ,2)  ~2+state (i ,3)  ~2)*  .  .  . 
(target (i , 1) ~2+target (i , 2) ~2+target (i ,3) ~2) ) ; 
num=state(i , 1) ^target (i , l)+state (i , 2) ^target (i ,2)+. . . 
state(i,3)*target(i,3) ; 

r_target_2= (target (i , 1) "2+target (i , 2) ~2+target (i ,3) ~2) ; 
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H=  [-stated ,  1)  *num*r_target_2/denunT3+target  (i ,  1)  /denum 
-state (i , 2) *num*r_target_2/denunT3+target (i , 2) / denum 
-state (i , 3) *num*r_target_2/denunT3+target (i , 3) / denum 

0  .  .  . 

0  .  .  . 

0]; 


°/„  computing  residuals 

G(i , 1 : l)=state (i , 1 : 3) ^target (i , 1 : 3) V (norm (state (i ,1:3))*. . . 

norm (target (i , 1:3))); 

rsdl (i , 1 : l)=z_i (i , clmn) — G(i ,1:1) ; 


T=H*phi ; 


P_inv_l=T) *inv(Q (i , clmn) ) *T ; 
rsdl_sum_l=T’ *inv(Q (i , clmn) ) *rsdl (i , 1:1); 


if  i>lastdata 
P_inv_l=0; 
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rsdl_sum_l=0 ; 

end 

P_inv=P_inv+P_inv_l ; 
rsdl_sum=  rsdl_sum+rsdl_sum_l ; 

end 

P=inv(P_inv) ; 
delta_x=P*rsdl_sum; 
dbstop  if  warning 
state_new=state (1,1:6) ’ +delta_x ; 
clear  y  rsdl  state 

end 

[t , x] =ode45 (@state_eom , t , state_new) ; 

state_last ( : , 1 : 6)=x( : , 1 : 6) ; 
clear  x 
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A. 2.3  The  Function  to  Compute  the  Estimated  State  at  the  Impact  Time  with 


Respect  to  the  Body  Frame. 

function  [pos ,pos_cross ,pos_radial] =compute_r (state_last , clmn, target) 

du=6378. 135; 

tu=sqrt (6378 . 135~3/3 . 98601e5) ; 
mu=l ; 


bl=target (length ( st at e_last) ,4 : 6) /norm (target ( length (state_last) ,4:6)) ; 
b3=cross (bl , target (length(state_last) ,1:3)... 

/norm(target (length(state_last) ,1:3))); 
b3=b3/norm(b3) ; 
b2=cross (bl ,b3) ; 

R_ib=  [bl  ’  b2’  b3 5 ] ; 

R_bi=R_ib’ ; 

state_body=R_bi* (state_last (length(state_last) , 1 : 3) - . . . 
target (length(state_last) , 1 : 3) ) J ; 

pos=norm(state_body(2) , state_body (3) ) ; 
pos_cross=state_body (2) ; 
pos_radial=state_body (3) ; 
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figured) 

line( [0  0] ,  [0  0] , ’Marker’ , ’d’ , ’Markerf acecolor ’ , ’r’ , ’LineStyle’ . . . 
, ’ none ’ , ’ color ’ , ’ r ’ ) 

line ( [state_body (2) *du  0] , [state_body (3) *du  0] , ’Marker’ , .. . 

’ LineStyle ’ , ’ none ’ , ’ color ’ , ’ b ’ ) 

clear  state_last 

A. 2. 4  The  Function  to  Iterate  the  State  and  the  $  Matrix  in  Time. 

function  y_dot=y_eom_l (t ,y) 
mu=l ; 


X=y(37); 

Y=y(38) ; 

Z=y(39); 

r=sqrt(X~2+Y~2+Z~2) ; 

B=[zeros(3)  eye (3) ; (-mu/r“3) *eye (3)  zeros(3)] ; 


A_rr=[(-mu/ (r~3))+((3*mu*X~2)/ (r"5))  (3*mu*X*Y)/ (r~5)  (3*mu*X*Z)/ (r~5) ; . . . 
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(3*mu*X*Y)/ (r~5)  (-mu/(r~3) )+( (3*mu*Y~2)/ (r~5))  (3*mu*Y*Z)/ (r~5) ; . . . 
(3*mu*X*Z)/ (r~5)  (3*mu*Y*Z)/ (r~5)  (-mu/ (r~3) )+( (3*mu*Z~2)/ (r~5))] ; 

A= [zeros (3)  eye(3);A_rr  zeros (3) ] ; 

AA=  [A ( 1 , 1) *eye(6)  A ( 1 , 2) *eye (6)  A ( 1 , 3) *eye (6)  A(l,4)*eye(6) . . . 

A(l,5)*eye(6)  A ( 1 , 6) *eye (6) ; . . . 

A(2 , 1) *eye (6)  A(2 , 2) *eye (6)  A(2,3)*eye(6)  A(2,4)*eye(6) . . . 

A(2 , 5) *eye (6)  A(2 , 6) *eye (6) ; . . . 

A(3, l)*eye(6)  A(3 , 2) *eye (6)  A(3,3)*eye(6)  A(3,4)*eye(6) . . . 

A(3,5)*eye(6)  A(3 , 6) *eye (6) ; . . . 

A(4, l)*eye(6)  A(4,2)*eye(6)  A(4,3)*eye(6)  A(4,4)*eye(6) . . . 

A(4,5)*eye(6)  A(4,6)*eye(6) ; . . . 

A(5 , 1) *eye (6)  A(5 , 2) *eye (6)  A(5,3)*eye(6)  A(5,4)*eye(6) . . . 

A(5 , 5) *eye (6)  A(5 , 6) *eye (6) ; . . . 

A(6 , 1) *eye (6)  A(6 , 2) *eye (6)  A(6,3)*eye(6)  A(6,4)*eye(6) . . . 

A(6 , 5) *eye (6)  A(6 , 6) *eye (6) ; ] ; 

AAA= [AA  zeros (36, 6) ; zeros (6, 36)  B] ; 

y_dot=AAA*y ; 

end 
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A. 2. 5  The  Function  to  Compute  the  Mean  and  the  Standard  Deviation  Values. 


function  [mean, stdev]  =  stat(x) 
n  =  length (x); 
mean  =  sum(x)/n; 

stdev  =  sqrt (sum( (x-mean) . ~2/n) ) ; 
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